module GOLD_initialization
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, April 1994 - June 2002                         *
!*                                                                     *
!*    This subroutine initializes the fields for the simulations.      *
!*  The one argument passed to initialize, Time, is set to the         *
!*  current time of the simulation.  The fields which are initialized  *
!*  here are:                                                          *
!*    u - Zonal velocity in m s-1.                                     *
!*    v - Meridional velocity in m s-1.                                *
!*    h - Layer thickness in m.  (Must be positive.)                   *
!*    D - Basin depth in m.  (Must be positive.)                       *
!*    f - The Coriolis parameter, in s-1.                              *
!*    g - The reduced gravity at each interface, in m s-2.             *
!*    Rlay - Layer potential density (coordinate variable) in kg m-3.  *
!*  If ENABLE_THERMODYNAMICS is defined:                               *
!*    T - Temperature in C.                                            *
!*    S - Salinity in psu.                                             *
!*  If SPONGE is defined:                                              *
!*    A series of subroutine calls are made to set up the damping      *
!*    rates and reference profiles for all variables that are damped   *
!*    in the sponge.                                                   *
!*  Any user provided tracer code is also first linked through this    *
!*  subroutine.                                                        *
!*                                                                     *
!*    Forcing-related fields (taux, tauy, buoy, ustar, etc.) are set   *
!*  in GOLD_surface_forcing.F90.                                       *
!*                                                                     *
!*    These variables are all set in the set of subroutines (in this   *
!*  file): initialize_topography, initialize_thickness,                *
!*  initialize_velocity, initialize_temp_sal, initialize_sponges, and  *
!*  set_coordinate.                                                    *
!*                                                                     *
!*  Macros written all in capital letters are defined in GOLD_memory.h.*
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, f                                     *
!*    j+1  > o > o >   At ^:  v, tauy                                  *
!*    j    x ^ x ^ x   At >:  u, taux                                  *
!*    j    > o > o >   At o:  h, D, buoy, tr, T, S, Rml, ustar         *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use GOLD_checksums, only : hchksum, qchksum, uchksum, vchksum, chksum
use GOLD_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast
use GOLD_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE
use GOLD_error_handler, only : GOLD_mesg, GOLD_error, FATAL, WARNING, is_root_pe
use GOLD_file_parser, only : read_param, get_param, open_param_file, param_file_type
use GOLD_file_parser, only : uppercase, log_param, log_version
use GOLD_grid, only : ocean_grid_type
use GOLD_interface_heights, only : find_eta
use GOLD_io, only : close_file, create_file, fieldtype, file_exists
use GOLD_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE
use GOLD_io, only : slasher, vardesc, write_field
use GOLD_io, only : EAST_FACE, NORTH_FACE
use GOLD_grid_initialize, only : initialize_masks, set_grid_metrics
use GOLD_restart, only : restore_state, GOLD_restart_CS
use GOLD_sponge, only : set_up_sponge_field, set_up_sponge_ML_density
use GOLD_sponge, only : initialize_sponge, sponge_CS
use GOLD_time_manager, only : time_type, set_time
use GOLD_tracer, only : add_tracer_OBC_values, advect_tracer_CS
use GOLD_variables, only : thermo_var_ptrs, directories, ocean_OBC_type
use GOLD_variables, only : OBC_NONE, OBC_SIMPLE, OBC_FLATHER_E, OBC_FLATHER_W
use GOLD_variables, only : OBC_FLATHER_N, OBC_FLATHER_S
use GOLD_EOS, only : calculate_density, calculate_density_derivs, EOS_type
use GOLD_EOS, only : int_specific_vol_dp
use user_initialization, only : user_set_coord, user_initialize_topography
use user_initialization, only : user_initialize_thickness, user_initialize_velocity
use user_initialization, only : user_init_temperature_salinity
use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges
use DOME_initialization, only : DOME_initialize_thickness
use DOME_initialization, only : DOME_initialize_topography
use DOME_initialization, only : DOME_set_Open_Bdry_Conds
use DOME_initialization, only : DOME_initialize_sponges
use benchmark_initialization, only : benchmark_initialize_thickness
use benchmark_initialization, only : benchmark_initialize_topography
use benchmark_initialization, only : benchmark_init_temperature_salinity
use circle_obcs_initialization, only : circle_obcs_initialize_thickness
use lock_exchange_initialization, only : lock_exchange_initialize_thickness
use external_gwave_initialization, only : external_gwave_initialize_thickness
use midas_vertmap, only : fill_miss_2d, find_interfaces, tracer_Z_init, meshgrid
use midas_vertmap, only : determine_temperature

use mpp_domains_mod, only : mpp_global_field, mpp_get_compute_domain
use horiz_interp_mod, only : horiz_interp_new, horiz_interp,horiz_interp_type
use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del

use netcdf

implicit none ; private

#include <GOLD_memory.h>

public GOLD_initialize, Get_GOLD_Input
public GOLD_initialize_topography, GOLD_initialize_rotation

! This structure is to simplify communication with the calling code.
type, public :: GOLD_initialization_struct
  type(advect_tracer_CS), pointer :: advect_tracer_CSp => NULL()
  type(sponge_CS), pointer :: sponge_CSp => NULL()
  type(ocean_OBC_type), pointer :: OBC => NULL()
end type GOLD_initialization_struct

contains

! -----------------------------------------------------------------------------
subroutine GOLD_initialize(u, v, h, tv, Time, G, PF, dirs, &
                              restart_CS, CS, Time_in)
  real, dimension(NXMEMQ_,NYMEM_,NZ_), intent(out)   :: u
  real, dimension(NXMEM_,NYMEMQ_,NZ_), intent(out)   :: v
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(out)   :: h
  type(thermo_var_ptrs),               intent(inout) :: tv
  type(time_type),                     intent(inout) :: Time
  type(ocean_grid_type),               intent(inout) :: G
  type(param_file_type),               intent(in)    :: PF
  type(directories),                   intent(in)    :: dirs
  type(GOLD_restart_CS),                pointer       :: restart_CS
  type(GOLD_initialization_struct),     intent(inout) :: CS
  type(time_type), optional,           intent(in) :: Time_in
! Arguments: u  - Zonal velocity, in m s-1.
!  (out)     v  - Meridional velocity, in m s-1.
!  (out)     h  - Layer thickness, in m.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (out)     Time    - Time at the start of the run segment.
!  (inout)   G       - The ocean's grid structure.
!  (in)      PF      - A structure indicating the open file to parse for
!                      model parameter values.
!  (in)      dirs    - A structure containing several relevant directory paths.
!  (inout)   restart_CS - A pointer to the restart control structure.
!  (inout)   CS      - A structure of pointers to be exchanged with GOLD.F90.
!  (in)      Time_in - Time at the start of the run segment. Time_in overrides
!                      any value set for Time.

  character(len=200) :: filename   ! The name of an input file.
  character(len=200) :: filename2  ! The name of an input files.
  character(len = 200) :: inputdir ! The directory where NetCDF input files are.
  character(len=200) :: config
  logical :: from_Z_file
  logical :: new_sim, write_geom
  logical :: use_temperature, use_sponge
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  logical :: depress_sfc ! If true, remove the mass that would be displaced
                         ! by a large surface pressure, such as with an ice
                         ! sheet.
  logical :: Analytic_FV_PGF, obsol_test
  logical :: apply_OBC_u, apply_OBC_v
  logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west
  logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south
  logical :: convert
  type(EOS_type), pointer :: eos => NULL()
  logical :: debug    ! indicates whether to write debugging output
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod = "GOLD_initialization" ! This module's name.
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  call GOLD_mesg(" GOLD_initialization.F90, GOLD_initialize: subroutine entered", 3)
  call log_version(PF, mod, version, tagname)

  new_sim = .false.
  if ((dirs%input_filename(1:1) == 'n') .and. &
      (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true.

  call get_param(PF, mod, "INPUTDIR", inputdir, &
         "The directory in which input files are found.", default=".")
  inputdir = slasher(inputdir)

  use_temperature = ASSOCIATED(tv%T)
  use_EOS = associated(tv%eqn_of_state)
  if (use_EOS) eos => tv%eqn_of_state

  call get_param(PF, mod, "DEBUG", debug, default=.false.)

! ====================================================================
!    Initialize fields that are time invariant - metrics, topography,
!  masks, vertical coordinate, Coriolis parameter.
! ====================================================================

! Set-up the layer densities, G%Rlay, and reduced gravities, G%g_prime.
  call get_param(PF, mod, "COORD_CONFIG", config, &
                 "This specifies how layers are to be defined: \n"//&
                 " \t file - read coordinate information from the file \n"//&
                 " \t\t specified by (COORD_FILE).\n"//&
                 " \t ts_ref - use reference temperature and salinity \n"//&
                 " \t ts_range - use range of temperature and salinity \n"//&
                 " \t\t (T_REF and S_REF) to determine surface density \n"//&
                 " \t\t and GINT calculate internal densities. \n"//&
                 " \t gprime - use reference density (RHO_0) for surface \n"//&
                 " \t\t density and GINT calculate internal densities. \n"//&
                 " \t ts_profile - use temperature and salinity profiles \n"//&
                 " \t\t (read from COORD_FILE) to set layer densities. \n"//&
                 " \t USER - call a user modified routine.", &
                 fail_if_missing=.true.)
  select case ( trim(config) )
    case ("gprime")
      call set_coord_from_gprime(G%Rlay, G%g_prime, G, PF)
    case ("layer_ref")
      call set_coord_from_layer_density(G%Rlay, G%g_prime, G, PF)
    case ("ts_ref")
      call set_coord_from_ts_ref(G%Rlay, G%g_prime, G, PF, eos, tv%P_Ref)
    case ("ts_profile")
      call set_coord_from_TS_profile(G%Rlay, G%g_prime, G, PF, eos, tv%P_Ref)
    case ("ts_range")
      call set_coord_from_TS_range(G%Rlay, G%g_prime, G, PF, eos, tv%P_Ref)

    case ("file")
      call set_coord_from_file(G%Rlay, G%g_prime, G, PF)
    case ("USER")
      call user_set_coord(G%Rlay, G%g_prime, G, PF, eos)
    case default ; call GOLD_error(FATAL,"GOLD_initialize: "// &
      "Unrecognized coordinate setup"//trim(config))
  end select
  if (debug) call chksum(G%Rlay, "GOLD_initialize: Rlay ", 1, nz)
  if (debug) call chksum(G%g_prime, "GOLD_initialize: g_prime ", 1, nz)

! Set up the parameters of the physical domain (i.e. the grid), G
  call set_grid_metrics(G, PF)

! Set up the bottom depth, G%D either analytically or from file
  call GOLD_initialize_topography(G%D, G, PF)

!    This call sets seamasks that prohibit flow over any point with  !
!  a bottom that is shallower than min_depth from PF.                !
  call initialize_masks(G, PF)
  if (debug) then
    call hchksum(G%D, 'GOLD_initialize: depth ', G, haloshift=1)
    call hchksum(G%hmask, 'GOLD_initialize: hmask ', G)
    call uchksum(G%umask, 'GOLD_initialize: umask ', G)
    call vchksum(G%vmask, 'GOLD_initialize: vmask ', G)
    call qchksum(G%qmask, 'GOLD_initialize: qmask ', G)
  endif

! Modulate geometric scales according to geography.
  call get_param(PF, mod, "CHANNEL_CONFIG", config, &
                 "A parameter that determines which set of channels are \n"//&
                 "restricted to specific  widths.  Options are:\n"//&
                 " \t none - All channels have the grid width.\n"//&
                 " \t global_1deg - Sets 16 specific channels appropriate \n"//&
                 " \t\t for a 1-degree model, as used in CM2G.\n"//&
                 " \t list - Read the channel locations and widths from a \n"//&
                 " \t\t text file, like GOLD_channel_list in the GOLD_SIS \n"//&
                 " \t\t test case.\n"//&
                 " \t file - Read open face widths everywhere from a \n"//&
                 " \t\t NetCDF file on the model grid.", &
                 default="none")
  select case ( trim(config) )
    case ("none")
    case ("list") ; call reset_face_lengths_list(G, PF)
    case ("file") ; call reset_face_lengths_file(G, PF)
    case ("global_1deg") ; call reset_face_lengths_named(G, PF, trim(config))
    case default ; call GOLD_error(FATAL, "GOLD_initialize: "// &
      "Unrecognized channel configuration "//trim(config))
  end select

!   This call sets the topography at velocity points.
  if (G%bathymetry_at_vel) then
    call get_param(PF, mod, "VELOCITY_DEPTH_CONFIG", config, &
                   "A string that determines how the topography is set at \n"//&
                   "velocity points. This may be 'min' or 'max'.", &
                   default="max")
    select case ( trim(config) )
      case ("max") ; call set_velocity_depth_max(G)
      case ("min") ; call set_velocity_depth_min(G)
      case default ; call GOLD_error(FATAL, "GOLD_initialize: "// &
        "Unrecognized velocity depth configuration "//trim(config))
    end select
  endif

!    Calculate the value of the Coriolis parameter at the latitude   !
!  of the q grid points, in s-1.
  call GOLD_initialize_rotation(G%f, G, PF)
  if (debug) then
    call qchksum(G%f, "GOLD_initialize: f ", G)
  endif

! Write out all of the grid data used by this run.
  call get_param(PF, mod, "ALWAYS_WRITE_GEOM", write_geom, &
                 "If true, write the geometry and vertical grid files \n"//&
                 "every time the model is run.  Otherwise, only write \n"//&
                 "them for new runs.", default=.true.)
  if (write_geom .or. new_sim) then
    call write_ocean_geometry_file(G, PF, dirs%output_directory)
    call write_vertgrid_file(G, PF, dirs%output_directory)
  endif

!====================================================================
!    Initialize temporally evolving fields, either as initial
!  conditions or by reading them from a restart (or saves) file.
!====================================================================

  if (new_sim) then
!  This block initializes all of the fields internally.              !
    call GOLD_mesg("Run initialized internally.", 3)

    if (present(Time_in)) Time = Time_in
    ! Otherwise leave Time at its input value.

    call get_param(PF, mod, "INIT_LAYERS_FROM_Z_FILE", from_Z_file, &
               "If true, intialize the layer thicknesses, temperatures, \n"//&
               "and salnities from a Z-space file on a latitude- \n"//&
               "longitude grid.", default=.false.)

    if (from_Z_file) then
!     Initialize thickness and T/S from z-coordinate data in a file.
      if (.NOT.use_temperature) call GOLD_error(FATAL,"GOLD_initialize : "//&
         "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")

      call GOLD_temp_salt_initialize_from_Z(h, tv, G, PF, dirs)
      call pass_var(h, G%Domain)
    else
!     Initialize thickness, h.
      call get_param(PF, mod, "THICKNESS_CONFIG", config, &
               "A string that determines how the initial layer \n"//&
               "thicknesses are specified for a new run: \n"//&
               " \t file - read interface heights from the file specified \n"//&
               " \t thickness_file - read thicknesses from the file specified \n"//&
               " \t\t by (THICKNESS_FILE).\n"//&
               " \t uniform - uniform thickness layers evenly distributed \n"//&
               " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
               " \t DOME - use a slope and channel configuration for the \n"//&
               " \t\t DOME sill-overflow test case. \n"//&
               " \t benchmark - use the benchmark test case thicknesses. \n"//&
               " \t search - search a density profile for the interface \n"//&
               " \t\t densities. This is not yet implemented. \n"//&
               " \t circle_obcs - the circle_obcs test case is used. \n"//&
               " \t USER - call a user modified routine.", &
               fail_if_missing=.true.)
      select case (trim(config))
         case ("file");    call initialize_thickness_from_file(h, G, PF, .false.)
         case ("thickness_file"); call initialize_thickness_from_file(h, G, PF, .true.)
         case ("uniform"); call initialize_thickness_uniform(h, G, PF)
         case ("DOME");    call DOME_initialize_thickness(h, G, PF)
         case ("benchmark"); call benchmark_initialize_thickness(h, G, PF, &
                                 tv%eqn_of_state, tv%P_Ref)
         case ("search");  call initialize_thickness_search
         case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, PF)
         case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, PF)
         case ("external_gwave"); call external_gwave_initialize_thickness(h, G, PF)
         case ("USER");    call user_initialize_thickness(h, G, PF,tv%T)
         case default ; call GOLD_error(FATAL,  "GOLD_initialize: "//&
              "Unrecognized layer thickness configuration "//trim(config))
      end select
      call pass_var(h, G%Domain)

!     Initialize temperature and salinity (T and S).
      if ( use_temperature ) then
        call get_param(PF, mod, "TS_CONFIG", config, &
               "A string that determines how the initial tempertures \n"//&
               "and salinities are specified for a new run: \n"//&
               " \t file - read velocities from the file specified \n"//&
               " \t\t by (TS_FILE). \n"//&
               " \t fit - find the temperatures that are consistent with \n"//&
               " \t\t the layer densities and salinity S_REF. \n"//&
               " \t TS_profile - use temperature and salinity profiles \n"//&
               " \t\t (read from TS_FILE) to set layer densities. \n"//&
               " \t benchmark - use the benchmark test case T & S. \n"//&
               " \t USER - call a user modified routine.", &
               fail_if_missing=.true.)
        select case (trim(config))
          case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, G, PF, eos, tv%P_Ref)
          case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, G, PF)
          case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
                                       G, PF, eos, tv%P_Ref)
          case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, G, PF)
          case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, G, PF, eos)
          case default ; call GOLD_error(FATAL,  "GOLD_initialize: "//&
                 "Unrecognized Temp & salt configuration "//trim(config))
        end select
      endif
    endif  ! not from_Z_file.

    if (debug) then
      call hchksum(h, "GOLD_initialize: h ", G)
      if ( use_temperature ) call hchksum(tv%T, "GOLD_initialize: T ", G)
      if ( use_temperature ) call hchksum(tv%S, "GOLD_initialize: S ", G)
    endif

!   Initialize velocity components, u and v
    call get_param(PF, mod, "VELOCITY_CONFIG", config, &
         "A string that determines how the initial velocities \n"//&
         "are specified for a new run: \n"//&
         " \t file - read velocities from the file specified \n"//&
         " \t\t by (VELOCITY_FILE). \n"//&
         " \t zero - the fluid is initially at rest. \n"//&
         " \t uniform - the flow is uniform (determined by\n"//&
         " \t\t paremters TORUS_U and TORUS_V).\n"//&
         " \t USER - call a user modified routine.", default="zero")
    select case (trim(config))
       case ("file"); call initialize_velocity_from_file(u, v, G, PF)
       case ("zero"); call initialize_velocity_zero(u, v, G, PF)
       case ("uniform"); call initialize_velocity_uniform(u, v, G, PF)
       case ("circular"); call initialize_velocity_circular(u, v, G, PF)
       case ("USER"); call user_initialize_velocity(u, v, G, PF)
       case default ; call GOLD_error(FATAL,  "GOLD_initialize: "//&
            "Unrecognized velocity configuration "//trim(config))
    end select

    call pass_vector(u, v, G%Domain)
    if (debug) call uchksum(u, "GOLD_initialize: u ", G)
    if (debug) call vchksum(v, "GOLD_initialize: v ", G)




!   Optionally convert the thicknesses from m to kg m-2.  This is particularly
! useful in a non-Boussinesq model.
    call get_param(PF, mod, "CONVERT_THICKNESS_UNITS", convert, &
                 "If true,  convert the thickness initial conditions from \n"//&
                 "units of m to kg m-2 or vice versa, depending on whether \n"//&
                 "BOUSSINESQ is defined. This does not apply if a restart \n"//&
                 "file is read.", default=.false.)
    if (convert) call convert_thickness(h, G, PF, tv)

!  Remove the mass that would be displaced by an ice shelf or inverse barometer.
    call get_param(PF, mod, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
                 "If true,  depress the initial surface to avoid huge \n"//&
                 "tsunamis when a large surface pressure is applied.", &
                 default=.false.)
    if (depress_sfc) call depress_surface(h, G, PF, tv)

  else ! Previous block for new_sim=.T., this block restores state
!    This line calls a subroutine that reads the initial conditions  !
!  from a previously generated file.                                 !
    call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, &
                       G, restart_CS)
    if (present(Time_in)) Time = Time_in
  endif

  call get_param(PF, mod, "SPONGE", use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. \n"//&
                 "The exact location and properties of those sponges are \n"//&
                 "specified via SPONGE_CONFIG.", default=.false.)
  if ( use_sponge ) then
! The 3 arguments here are (1) a flag indicating whether the sponge  !
! values are to be read from files, (2) the name of a file containing!
! the state toward which the model is damped, and (3) the file in    !
! which the 2-D damping rate field can be found.                     !
    call get_param(PF, mod, "SPONGE_CONFIG", config, &
                 "A string that sets how the sponges are configured: \n"//&
                 " \t file - read sponge properties from the file \n"//&
                 " \t\t specified by (SPONGE_FILE).\n"//&
                 " \t DOME - use a slope and channel configuration for the \n"//&
                 " \t\t DOME sill-overflow test case. \n"//&
                 " \t USER - call a user modified routine.", default="file")

    select case (trim(config))
      case ("DOME"); call DOME_initialize_sponges(G, tv, PF, CS%sponge_CSp)
      case ("USER"); call user_initialize_sponges(G, use_temperature, tv, &
                                                  PF, CS%sponge_CSp, h)
      case ("file"); call initialize_sponges_file(G, use_temperature, tv, &
                                                  PF, CS%sponge_CSp)
      case default ; call GOLD_error(FATAL,  "GOLD_initialize: "//&
             "Unrecognized sponge configuration "//trim(config))
    end select

  endif

! This subroutine call sets optional open boundary conditions.
  call get_param(PF, mod, "APPLY_OBC_U", apply_OBC_u, &
                 "If true, open boundary conditions may be set at some \n"//&
                 "u-points, with the configuration controlled by OBC_CONFIG", &
                 default=.false.)
  call get_param(PF, mod, "APPLY_OBC_V", apply_OBC_v, &
                 "If true, open boundary conditions may be set at some \n"//&
                 "v-points, with the configuration controlled by OBC_CONFIG", &
                 default=.false.)
  if (apply_OBC_u .or. apply_OBC_v) then 
    call get_param(PF, mod, "OBC_CONFIG", config, &
                 "A string that sets how the open boundary conditions are \n"//&
                 " configured: \n"//&
                 " \t DOME - use a slope and channel configuration for the \n"//&
                 " \t\t DOME sill-overflow test case. \n"//&
                 " \t USER - call a user modified routine.", default="file", &
                 fail_if_missing=.true.)
    if (trim(config) == "DOME") then
      call DOME_set_Open_Bdry_Conds(CS%OBC, tv, G, PF, CS%advect_tracer_CSp)
    elseif (trim(config) == "USER") then
      call user_set_Open_Bdry_Conds(CS%OBC, tv, G, PF, CS%advect_tracer_CSp)
    else
      call GOLD_error(FATAL, "The open boundary conditions specified by "//&
              "OBC_CONFIG = "//trim(config)//" have not been fully implemented.")
      call set_Open_Bdry_Conds(CS%OBC, tv, G, PF, CS%advect_tracer_CSp)
    endif
  endif

  call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,&
                 "Apply a Flather open boundary condition on the eastern \n"//&
                 "side of the global domain", default=.false.)
  call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,&
                 "Apply a Flather open boundary condition on the western \n"//&
                 "side of the global domain", default=.false.)
  call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,&
                 "Apply a Flather open boundary condition on the northern \n"//&
                 "side of the global domain", default=.false.)
  call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,&
                 "Apply a Flather open boundary condition on the southern \n"//&
                 "side of the global domain", default=.false.)
  if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) then
    call set_Flather_Bdry_Conds(CS%OBC, tv, h, G, PF, CS%advect_tracer_CSp)
  endif

  call GOLD_mesg(" GOLD_initialization.F90, GOLD_initialize: complete", 3)

end subroutine GOLD_initialize
! -----------------------------------------------------------------------------


! -----------------------------------------------------------------------------
subroutine set_coord_from_gprime(Rlay, g_prime, G, param_file)
  real, dimension(:), intent(out)   :: Rlay, g_prime
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
! Arguments: Rlay - the layers' target coordinate values (potential density).
!  (out)     g_prime - the reduced gravity across the interfaces, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

! This subroutine sets the layer densities (Rlay) and the interface  !
! reduced gravities (g).                                             !
  real :: g_int   ! Reduced gravities across the internal interfaces, in m s-2.
  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
  character(len=40)  :: mod = "set_coord_from_gprime" ! This subroutine's name.
  integer :: k, nz
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, set_coord_from_gprime: setting coordinate", 5)

  call get_param(param_file, mod, "GFS" , g_fs, &
                 "The reduced gravity at the free surface.", units="m s-2", &
                 default=G%g_Earth)
  call get_param(param_file, mod, "GINT", g_int, &
                 "The reduced gravity across internal interfaces.", &
                 units="m s-2", fail_if_missing=.true.)

  g_prime(1) = g_fs
  do k=2,nz ; g_prime(k) = g_int ; enddo
  Rlay(1) = G%Rho0
  do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(G%Rho0/G%g_Earth) ; enddo


end subroutine set_coord_from_gprime
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_coord_from_layer_density(Rlay, g_prime, G, param_file)
  real, dimension(:), intent(out)   :: Rlay, g_prime
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
! Arguments: Rlay - the layers' target coordinate values (potential density).
!  (out)     g_prime - the reduced gravity across the interfaces, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

! This subroutine sets the layer densities (Rlay) and the interface  !
! reduced gravities (g).                                             !
  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
  real :: Rlay_Ref! The surface layer's target density, in kg m-3.
  real :: RLay_range ! The range of densities, in kg m-3.
  character(len=40)  :: mod = "set_coord_from_layer_density" ! This subroutine's name.
  integer :: k, nz
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, set_coord_from_layer_density: setting coordinate", 5)

  call get_param(param_file, mod, "GFS", g_fs, &
                 "The reduced gravity at the free surface.", units="m s-2", &
                 default=G%g_Earth)
  call get_param(param_file, mod, "LIGHTEST_DENSITY", Rlay_Ref, &
                 "The reference potential density used for layer 1.", &
                 units="kg m-3", default=G%Rho0)
  call get_param(param_file, mod, "DENSITY_RANGE", Rlay_range, &
                 "The range of reference potential densities in the layers.", &
                 units="kg m-3", default=2.0)

  g_prime(1) = g_fs
  Rlay(1) = Rlay_Ref
  do k=2,nz
     Rlay(k) = Rlay(k-1) + RLay_range/(real(nz-1))
  enddo
!    These statements set the interface reduced gravities.           !
  do k=2,nz
     g_prime(k) = (G%g_Earth/G%Rho0) * (Rlay(k) - Rlay(k-1))
  enddo

end subroutine set_coord_from_layer_density
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_coord_from_TS_ref(Rlay, g_prime, G, param_file, eqn_of_state, &
                                 P_Ref)
  real, dimension(:),    intent(out) :: Rlay, g_prime
  type(ocean_grid_type), intent(in)  :: G
  type(param_file_type), intent(in)  :: param_file
  type(EOS_type),        pointer     :: eqn_of_state
  real,                  intent(in)  :: P_Ref
! Arguments: Rlay - the layers' target coordinate values (potential density).
!  (out)     g_prime - the reduced gravity across the interfaces, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      eqn_of_state - integer selecting the equation of state
!  (in)      P_Ref - The coordinate-density reference pressure in Pa.

! This subroutine sets the layer densities (Rlay) and the interface  !
! reduced gravities (g).                                             !
  real :: T_ref   ! Reference temperature
  real :: S_ref   ! Reference salinity
  real :: g_int   ! Reduced gravities across the internal interfaces, in m s-2.
  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
  character(len=40)  :: mod = "set_coord_from_TS_ref" ! This subroutine's name.
  integer :: k, nz
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, set_coord_from_TS_ref: setting coordinate", 5)

  call get_param(param_file, mod, "T_REF", T_Ref, &
                 "The initial temperature of the lightest layer.", units="degC", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "S_REF", S_Ref, &
                 "The initial salinities.", units="PSU", default=35.0)
  call get_param(param_file, mod, "GFS", g_fs, &
                 "The reduced gravity at the free surface.", units="m s-2", &
                 default=G%g_Earth)
  call get_param(param_file, mod, "GINT", g_int, &
                 "The reduced gravity across internal interfaces.", &
                 units="m s-2", fail_if_missing=.true.)
                                      !
!    These statements set the interface reduced gravities.           !
  g_prime(1) = g_fs
  do k=2,nz ; g_prime(k) = g_int ; enddo

!    The uppermost layer's density is set here.  Subsequent layers'  !
!  densities are determined from this value and the g values.        !
!        T0 = 28.228 ; S0 = 34.5848 ; Pref = P_Ref
  call calculate_density(T_ref, S_ref, P_ref, Rlay(1), 1,1, eqn_of_state)

!    These statements set the layer densities.                       !
  do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(G%Rho0/G%g_Earth) ; enddo


end subroutine set_coord_from_TS_ref
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_coord_from_TS_profile(Rlay, g_prime, G, param_file, &
                                     eqn_of_state, P_Ref)
  real, dimension(:),    intent(out) :: Rlay, g_prime
  type(ocean_grid_type), intent(in)  :: G
  type(param_file_type), intent(in)  :: param_file
  type(EOS_type),        pointer     :: eqn_of_state
  real,                  intent(in)  :: P_Ref
! Arguments: Rlay - the layers' target coordinate values (potential density).
!  (out)     g_prime - the reduced gravity across the interfaces, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      eqn_of_state - integer that selects equation of state
!  (in)      P_Ref - The coordinate-density reference pressure in Pa.

! This subroutine sets the layer densities (Rlay) and the interface  !
! reduced gravities (g).                                             !
  real, dimension(SZK_(G)) :: T0, S0,  Pref
  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
  integer :: k, nz
  character(len=40)  :: mod = "set_coord_from_TS_profile" ! This subroutine's name.
  character(len=200) :: filename, coord_file, inputdir ! Strings for file/path
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, set_coord_from_TS_profile: setting coordinate", 5)

  call get_param(param_file, mod, "GFS", g_fs, &
                 "The reduced gravity at the free surface.", units="m s-2", &
                 default=G%g_Earth)
  call get_param(param_file, mod, "COORD_FILE", coord_file, &
                 "The file from which the coordinate temperatures and \n"//&
                 "salnities are read.", fail_if_missing=.true.)

  call get_param(param_file,  mod, "INPUTDIR", inputdir, default=".")
  filename = trim(slasher(inputdir))//trim(coord_file)
  call log_param(param_file, mod, "INPUTDIR/COORD_FILE", filename)

  call read_data(filename,"PTEMP",T0(:),domain=G%Domain%mpp_domain)
  call read_data(filename,"SALT",S0(:),domain=G%Domain%mpp_domain)

  if (.not.file_exists(filename)) call GOLD_error(FATAL, &
      " set_coord_from_TS_profile: Unable to open " //trim(filename))
!    These statements set the interface reduced gravities.           !
  g_prime(1) = g_fs
  do k=1,nz ; Pref(k) = P_ref ; enddo
  call calculate_density(T0, S0, Pref, Rlay, 1,nz,eqn_of_state)
  do k=2,nz; g_prime(k) = (G%g_Earth/G%Rho0) * (Rlay(k) - Rlay(k-1)); enddo

end subroutine set_coord_from_TS_profile
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_coord_from_TS_range(Rlay, g_prime, G, param_file, &
                                     eqn_of_state, P_Ref)
  real, dimension(:),    intent(out) :: Rlay, g_prime
  type(ocean_grid_type), intent(in)  :: G
  type(param_file_type), intent(in)  :: param_file
  type(EOS_type),        pointer     :: eqn_of_state
  real,                  intent(in)  :: P_Ref
! Arguments: Rlay - the layers' target coordinate values (potential density).
!  (out)     g_prime - the reduced gravity across the interfaces, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      eqn_of_state - integer that selects equation of state
!  (in)      P_Ref - The coordinate-density reference pressure in Pa.

! This subroutine sets the layer densities (Rlay) and the interface  !
! reduced gravities (g).                                             !
  real, dimension(SZK_(G)) :: T0, S0,  Pref
  real :: S_Ref, S_Light, S_Dense ! Salnity range parameters in PSU.
  real :: T_Ref, T_Light, T_Dense ! Temperature range parameters in dec C.
  real :: res_rat ! The ratio of density space resolution in the denser part
                  ! of the range to that in the lighter part of the range.
                  ! Setting this greater than 1 increases the resolution for
                  ! the denser water.
  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
  real :: a1, frac_dense, k_frac
  integer :: k, nz, k_light
  character(len=40)  :: mod = "set_coord_from_TS_range" ! This subroutine's name.
  character(len=200) :: filename, coord_file, inputdir ! Strings for file/path
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, set_coord_from_TS_range: setting coordinate", 5)

  call get_param(param_file, mod, "T_REF", T_Ref, &
                 "The default initial temperatures.", units="degC", default=10.0)
  call get_param(param_file, mod, "TS_RANGE_T_LIGHT", T_Light, &
                 "The initial temperature of the lightest layer when \n"//&
                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref)
  call get_param(param_file, mod, "TS_RANGE_T_DENSE", T_Dense, &
                 "The initial temperature of the densest layer when \n"//&
                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref)

  call get_param(param_file, mod, "S_REF", S_Ref, &
                 "The default initial salinities.", units="PSU", default=35.0)
  call get_param(param_file, mod, "TS_RANGE_S_LIGHT", S_Light, &
                 "The initial lightest salinities when COORD_CONFIG \n"//&
                 "is set to ts_range.", default = S_Ref, units="PSU")
  call get_param(param_file, mod, "TS_RANGE_S_DENSE", S_Dense, &
                 "The initial densest salinities when COORD_CONFIG \n"//&
                 "is set to ts_range.", default = S_Ref, units="PSU")
		 
  call get_param(param_file, mod, "TS_RANGE_RESOLN_RATIO", res_rat, &
  		           "The ratio of density space resolution in the densest \n"//&
                 "part of the range to that in the lightest part of the \n"//&
                 "range when COORD_CONFIG is set to ts_range. Values \n"//&
                 "greater than 1 increase the resolution of the denser water.",&
                 default=1.0, units="nondim")

  call get_param(param_file, mod, "GFS", g_fs, &
                 "The reduced gravity at the free surface.", units="m s-2", &
                 default=G%g_Earth)

  k_light = G%nk_rho_varies + 1

  ! Set T0(k) to range from T_LIGHT to T_DENSE, and simliarly for S0(k).
  T0(k_light) = T_light ; S0(k_light) = S_light
  a1 = 2.0 * res_rat / (1.0 + res_rat)
  do k=k_light+1,nz
    k_frac = real(k-k_light)/real(nz-k_light)
    frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2
    T0(k) = frac_dense * (T_Dense - T_Light) + T_Light
    S0(k) = frac_dense * (S_Dense - S_Light) + S_Light
  enddo

  g_prime(1) = g_fs
  do k=1,nz ; Pref(k) = P_ref ; enddo
  call calculate_density(T0, S0, Pref, Rlay, k_light,nz-k_light+1,eqn_of_state)
  ! Extrapolate target densities for the variable density mixed and buffer layers.
  do k=k_light-1,1,-1
    Rlay(k) = 2.0*Rlay(k+1) -  Rlay(k+2)
  enddo
  do k=2,nz; g_prime(k) = (G%g_Earth/G%Rho0) * (Rlay(k) - Rlay(k-1)); enddo

end subroutine set_coord_from_TS_range
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_coord_from_file(Rlay, g_prime, G, param_file)
  real, dimension(:), intent(out)   :: Rlay, g_prime
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
! Arguments: Rlay - the layers' target coordinate values (potential density).
!  (out)     g_prime - the reduced gravity across the interfaces, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

! This subroutine sets the layer densities (Rlay) and the interface  !
! reduced gravities (g).                                             !
  real :: g_fs    ! Reduced gravity across the free surface, in m s-2.
  integer :: k, nz
  character(len=40)  :: mod = "set_coord_from_file" ! This subroutine's name.
  character(len=40)  :: coord_var
  character(len=200) :: filename,coord_file,inputdir ! Strings for file/path
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, set_coord_from_file: setting coordinate", 5)

  call get_param(param_file, mod, "GFS", g_fs, &
                 "The reduced gravity at the free surface.", units="m s-2", &
                 default=G%g_Earth)
  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mod, "COORD_FILE", coord_file, &
                 "The file from which the coordinate densities are read.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "COORD_VAR", coord_var, &
                 "The variable in COORD_FILE that is to be used for the \n"//&
                 "coordinate densities.", default="Layer")
  filename = trim(inputdir)//trim(coord_file)
  call log_param(param_file, mod, "INPUTDIR/COORD_FILE", filename)
  if (.not.file_exists(filename)) call GOLD_error(FATAL, &
      " set_coord_from_file: Unable to open "//trim(filename))

  call read_axis_data(filename, coord_var, Rlay)
  g_prime(1) = g_fs
  do k=2,nz ; g_prime(k) = (G%g_Earth/G%Rho0) * (Rlay(k) - Rlay(k-1)) ; enddo
  do k=1,nz ; if (g_prime(k) <= 0.0) then
    call GOLD_error(FATAL, "GOLD_initialization set_coord_from_file: "//&
       "Zero or negative g_primes read from variable "//"Layer"//" in file "//&
       trim(filename))
  endif ; enddo

end subroutine set_coord_from_file
! -----------------------------------------------------------------------------

subroutine GOLD_initialize_rotation(f, G, PF)
  type(ocean_grid_type),            intent(in)  :: G
  real, dimension(G%Isdq:G%Iedq,G%Jsdq:G%Jedq), intent(out) :: f
  type(param_file_type),            intent(in)  :: PF
! Arguments: f  - the Coriolis parameter in s-1. Intent out.
!  (in)      G  - The ocean's grid structure.
!  (in)      PF - A structure indicating the open file to parse for
!                         model parameter values.

!   This subroutine makes the appropriate call to set up the Coriolis parameter.
! This is a separate subroutine so that it can be made public and shared with
! the ice-sheet code or other components.
! Set up the Coriolis parameter, G%f, either analytically or from file
  character(len=40)  :: mod = "GOLD_initialize_rotation" ! This subroutine's name.
  character(len=200) :: config

  call get_param(PF, mod, "ROTATION", config, &
                 "This specifies how the Coriolis parameter is specified: \n"//&
                 " \t 2omegasinlat - Use twice the planetary rotation rate \n"//&
                 " \t\t times the sine of latitude.\n"//&
                 " \t betaplane - Use a beta-plane or f-plane. \n"//&
                 " \t USER - call a user modified routine.", &
                 default="2omegasinlat")
  select case (trim(config))
    case ("2omegasinlat"); call set_rotation_planetary(f, G, PF)
    case ("beta"); call set_rotation_beta_plane(f, G, PF)
    case ("betaplane"); call set_rotation_beta_plane(f, G, PF)
   !case ("nonrotating") ! Note from AJA: Missing case?
    case default ; call GOLD_error(FATAL,"GOLD_initialize: "// &
      "Unrecognized rotation setup "//trim(config))
  end select
end subroutine GOLD_initialize_rotation

subroutine GOLD_initialize_topography(D, G, PF)
  real, intent(out), dimension(NXMEM_,NYMEM_) :: D
  type(ocean_grid_type), intent(in)           :: G
  type(param_file_type), intent(in)           :: PF
! Arguments: D  - the bottom depth in m. Intent out.
!  (in)      G  - The ocean's grid structure.
!  (in)      PF - A structure indicating the open file to parse for
!                         model parameter values.

!  This subroutine makes the appropriate call to set up the bottom depth.
!  This is a separate subroutine so that it can be made public and shared with
!  the ice-sheet code or other components.
! Set up the bottom depth, G%D either analytically or from file
  character(len=40)  :: mod = "GOLD_initialize_topography" ! This subroutine's name.
  character(len=200) :: config

  call get_param(PF, mod, "TOPO_CONFIG", config, &
                 "This specifies how bathymetry is specified: \n"//&
                 " \t file - read bathymetric information from the file \n"//&
                 " \t\t specified by (TOPO_FILE).\n"//&
                 " \t flat - flat bottom set to MAXIMUM_DEPTH. \n"//&
                 " \t bowl - an analytically specified bowl-shaped basin \n"//&
                 " \t\t ranging between MAXIMUM_DEPTH and MINIMUM_DEPTH. \n"//&
                 " \t spoon - a similar shape to 'bowl', but with an vertical \n"//&
                 " \t\t wall at the southern face. \n"//&
                 " \t halfpipe - a zonally uniform channel with a half-sine \n"//&
                 " \t\t profile in the meridional direction. \n"//&
                 " \t DOME - use a slope and channel configuration for the \n"//&
                 " \t\t DOME sill-overflow test case. \n"//&
                 " \t benchmark - use the benchmark test case topography. \n"//&
                 " \t USER - call a user modified routine.", &
                 fail_if_missing=.true.)
  select case ( trim(config) )
    case ("file");      call initialize_topography_from_file(D, G, PF)
    case ("flat");      call initialize_topography_named(D, G, PF, config)
    case ("spoon");     call initialize_topography_named(D, G, PF, config)
    case ("bowl");      call initialize_topography_named(D, G, PF, config)
    case ("halfpipe");  call initialize_topography_named(D, G, PF, config)
    case ("DOME");      call DOME_initialize_topography(D, G, PF)
    case ("benchmark"); call benchmark_initialize_topography(D, G, PF)
    case ("USER");      call user_initialize_topography(D, G, PF)
    case default ;      call GOLD_error(FATAL,"GOLD_initialize_topography: "// &
      "Unrecognized topography setup '"//trim(config)//"'")
  end select
  if (trim(config) .ne. "DOME") then
    call limit_topography(D, G, PF)
  endif
  
end subroutine GOLD_initialize_topography

! -----------------------------------------------------------------------------
subroutine initialize_topography_from_file(D, G, param_file )
  real, intent(out), dimension(NXMEM_,NYMEM_) :: D
  type(ocean_grid_type), intent(in)           :: G
  type(param_file_type), intent(in)           :: param_file
! Arguments: D          - the bottom depth in m. Intent out.
!  (in)      G          - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!  This subroutine reads depths from a file and puts it into D(:,:) in m.
  character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
  character(len=200) :: topo_varname                  ! Variable name in file
  character(len=40)  :: mod = "initialize_topography_from_file" ! This subroutine's name.

  call GOLD_mesg("  GOLD_initialization.F90, initialize_topography_from_file: reading topography", 5)

  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mod, "TOPO_FILE", topo_file, &
                 "The file from which the bathymetry is read.", &
                 default="topog.nc")
  call get_param(param_file, mod, "TOPO_VARNAME", topo_varname, &
                 "The name of the bathymetry variable in TOPO_FILE.", &
                 default="depth")

  filename = trim(inputdir)//trim(topo_file)
  call log_param(param_file, mod, "INPUTDIR/TOPO_FILE", filename)

  if (.not.file_exists(filename, G%Domain)) call GOLD_error(FATAL, &
       " initialize_topography_from_file: Unable to open "//trim(filename))

  call read_data(filename,trim(topo_varname),D,domain=G%Domain%mpp_domain)

end subroutine initialize_topography_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_topography_named(D, G, param_file, topog_config)
  real, intent(out), dimension(NXMEM_,NYMEM_) :: D
  type(ocean_grid_type), intent(in)           :: G
  type(param_file_type), intent(in)           :: param_file
  character(len=*),      intent(in)           :: topog_config
! Arguments: D          - the bottom depth in m. Intent out.
!  (in)      G          - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      topog_config - The name of an idealized topographic configuration.

! This subroutine places the bottom depth in m into D(:,:), shaped in a spoon
  real :: min_depth, max_depth ! The minimum and maximum depths in m.
  real :: PI                   ! 3.1415926... calculated as 4*atan(1)
  real :: D0                   ! A constant to make the maximum     !
                               ! basin depth MAXIMUM_DEPTH.         !
  real :: expdecay             ! A decay scale of associated with   !
                               ! the sloping boundaries, in m.      !
  real :: Dedge                ! The depth in m at the basin edge.  !
  real :: southlat0, westlon0, lenlon0, lenlat0
  real :: south_lat, west_lon, len_lon, len_lat, Rad_earth
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, xhalo, yhalo
  character(len=40)  :: mod = "initialize_topography_named" ! This subroutine's name.
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed

  call GOLD_mesg("  GOLD_initialization.F90, initialize_topography_named: "//&
                 "setting topography "//trim(topog_config), 5)

  call get_param(param_file, mod, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0)
  call get_param(param_file, mod, "MAXIMUM_DEPTH", max_depth, &
                 "The maximum depth of the ocean.", units="m", &
                 fail_if_missing=.true.)

  ! These expressions force rounding of approximate values in a
  ! consistent way.
  xhalo = G%isc-G%isd ; yhalo = G%jsc-G%jsd
  southlat0 = ANInt(G%gridlatq(yhalo)*1024.0)/1024.0
  westlon0 = ANInt(G%gridlonq(xhalo)*1024.0)/1024.0
  lenlon0 = ANInt((G%gridlonq(G%Domain%nxtot+xhalo)-G%gridlonq(xhalo))*1024.0)/1024.0
  lenlat0 = ANInt((G%gridlatq(G%Domain%nytot+yhalo)-G%gridlatq(yhalo))*1024.0)/1024.0

  if (trim(topog_config) /= "flat") then
    call get_param(param_file, mod, "EDGE_DEPTH", Dedge, &
                   "The depth at the edge of one of the named topographies.", &
                   units="m", default=100.0)
    call get_param(param_file, mod, "SOUTHLAT", south_lat, &
                   "The southern latitude of the domain.", units="degrees", &
                   default=southlat0)
    call get_param(param_file, mod, "LENLAT", len_lat, &
                   "The latitudinal length of the domain.", units="degrees", &
                   default=lenlat0)
    call get_param(param_file, mod, "WESTLON", west_lon, &
                   "The western longitude of the domain.", units="degrees", &
                   default=westlon0)
    call get_param(param_file, mod, "LENLON", len_lon, &
                   "The longitudinal length of the domain.", units="degrees", &
                   default=lenlon0)
    call get_param(param_file, mod, "RAD_EARTH", Rad_Earth, &
                   "The radius of the Earth.", units="m", default=6.378e6)
    call get_param(param_file, mod, "TOPOG_SLOPE_SCALE", expdecay, &
                   "The exponential decay scale used in defining some of \n"//&
                   "the named topographies.", units="m", default=400000.0)
  endif


  PI = 4.0*atan(1.0)

  if (trim(topog_config) == "flat") then
    do i=is,ie ; do j=js,je ; D(i,j) = max_depth ; enddo ; enddo
  elseif (trim(topog_config) == "spoon") then
    D0 = (max_depth - Dedge) / &
             ((1.0 - exp(-0.5*len_lat*Rad_earth*PI/(180.0 *expdecay))) * &
              (1.0 - exp(-0.5*len_lat*Rad_earth*PI/(180.0 *expdecay))))
    do i=is,ie ; do j=js,je
  !  This sets a bowl shaped (sort of) bottom topography, with a       !
  !  maximum depth of max_depth.                                   !
      D(i,j) =  Dedge + D0 * &
             (sin(PI * (G%geolonh(i,j) - (west_lon)) / len_lon) * &
           (1.0 - exp((G%geolath(i,j) - (south_lat+len_lat))*Rad_earth*PI / &
                      (180.0*expdecay)) ))
    enddo ; enddo
  elseif (trim(topog_config) == "bowl") then
    D0 = (max_depth - Dedge) / &
             ((1.0 - exp(-0.5*len_lat*Rad_earth*PI/(180.0 *expdecay))) * &
              (1.0 - exp(-0.5*len_lat*Rad_earth*PI/(180.0 *expdecay))))

  !  This sets a bowl shaped (sort of) bottom topography, with a
  !  maximum depth of max_depth.
    do i=is,ie ; do j=js,je
      D(i,j) =  Dedge + D0 * &
             (sin(PI * (G%geolonh(i,j) - west_lon) / len_lon) * &
             ((1.0 - exp(-(G%geolath(i,j) - south_lat)*Rad_Earth*PI/ &
                          (180.0*expdecay))) * &
             (1.0 - exp((G%geolath(i,j) - (south_lat+len_lat))* &
                         Rad_Earth*PI/(180.0*expdecay)))))
    enddo ; enddo
  elseif (trim(topog_config) == "halfpipe") then
    D0 = max_depth - Dedge
    do i=is,ie ; do j=js,je
      D(i,j) =  Dedge + D0 * ABS(sin(PI*(G%geolath(i,j) - south_lat)/len_lat))
    enddo ; enddo
  else
    call GOLD_error(FATAL,"initialize_topography_named: "// &
      "Unrecognized topography name "//trim(topog_config))
  endif

  ! This is here just for safety.  Hopefully it doesn't do anything.
  do i=is,ie ; do j=js,je
    if (D(i,j) > max_depth) D(i,j) = max_depth
    if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth
  enddo ; enddo

end subroutine initialize_topography_named
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine limit_topography(D, G, param_file)
  real, intent(inout), dimension(NXMEM_,NYMEM_) :: D
  type(ocean_grid_type), intent(in)             :: G
  type(param_file_type), intent(in)             :: param_file
! Arguments: D          - the bottom depth in m. Intent in/out.
!  (in)      G          - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

! This subroutine ensures that    min_depth < D(x,y) < max_depth
  integer :: i, j
  character(len=40)  :: mod = "limit_topography" ! This subroutine's name.
  real :: min_depth, max_depth

  call GOLD_mesg("  GOLD_initialization.F90, limit_topography: limiting topography", 5)

  call get_param(param_file, mod, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0)
  call get_param(param_file, mod, "MAXIMUM_DEPTH", max_depth, &
                 "The maximum depth of the ocean.", units="m", default=1.0e6)

! Make sure that min_depth < D(x,y) < max_depth
  do j=G%jsd,G%jed ; do i=G%isd,G%ied
    D(i,j) = min( max( D(i,j), 0.5*min_depth ), max_depth )
  enddo ; enddo

end subroutine limit_topography
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_rotation_planetary(f, G, param_file)
  type(ocean_grid_type),            intent(in)  :: G
  real, dimension(G%Isdq:G%Iedq,G%Jsdq:G%Jedq), intent(out) :: f
  type(param_file_type),            intent(in)  :: param_file
! Arguments: f          - Coriolis parameter (vertical component) in s^-1
!     (in)   G          - grid type
!     (in)   param_file - parameter file type

! This subroutine sets up the Coriolis parameter for a sphere
  integer :: I, J
  real    :: PI, omega

  call GOLD_mesg("  GOLD_initialization.F90, set_rotation_planetary: setting f (Coriolis)", 5)

  call get_param(param_file, "set_rotation_planetary", "OMEGA", omega, &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5)
  PI = 4.0*atan(1.0)

  do I=G%Isdq,G%Iedq ; do J=G%Jsdq,G%Jedq
    f(I,J) = ( 2.0 * omega ) * sin( ( PI * G%geolatq(I,J) ) / 180.)
  enddo ; enddo

end subroutine set_rotation_planetary
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_rotation_beta_plane(f, G, param_file)
  type(ocean_grid_type),            intent(in)  :: G
  real, dimension(G%Isdq:G%Iedq,G%Jsdq:G%Jedq), intent(out) :: f
  type(param_file_type),            intent(in)  :: param_file
! Arguments: f          - Coriolis parameter (vertical component) in s^-1
!     (in)   G          - grid type
!     (in)   param_file - parameter file type

! This subroutine sets up the Coriolis parameter for a beta-plane
  integer :: I, J
  real    :: f_0, beta, y_scl, Rad_Earth, PI
  character(len=40)  :: mod = "set_rotation_beta_plane" ! This subroutine's name.
  character(len=200) :: axis_units

  call GOLD_mesg("  GOLD_initialization.F90, set_rotation_beta_plane: setting f (Coriolis)", 5)

  call get_param(param_file, mod, "F_0", f_0, &
                 "The reference value of the Coriolis parameter with the \n"//&
                 "betaplane option.", units="s-1", default=0.0)
  call get_param(param_file, mod, "BETA", beta, &
                 "The northward gradient of the Coriolis parameter with \n"//&
                 "the betaplane option.", units="m-1 s-1", default=0.0)
  call get_param(param_file, mod, "AXIS_UNITS", axis_units, default="degrees")

  PI = 4.0*atan(1.0)
  select case (axis_units(1:1))
    case ("d")
      call get_param(param_file, mod, "RAD_EARTH", Rad_Earth, &
                   "The radius of the Earth.", units="m", default=6.378e6)
      y_scl = Rad_Earth/PI
    case ("k"); y_scl = 1.E3
    case ("m"); y_scl = 1.
    case ("c"); y_scl = 1.E-2
    case default ; call GOLD_error(FATAL, &
      " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
  end select

  do I=G%Isdq,G%Iedq ; do J=G%Jsdq,G%Jedq
    f(I,J) = f_0 + beta * ( G%geolatq(I,J) * y_scl )
  enddo ; enddo

end subroutine set_rotation_beta_plane
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_thickness_from_file(h, G, param_file, file_has_thickness)
  real, intent(out), dimension(NXMEM_,NYMEM_, NZ_) :: h
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
  logical,               intent(in) :: file_has_thickness
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      file_has_thickness - If true, this file contains thicknesses;
!                                 otherwise it contains interface heights.

!  This subroutine reads the layer thicknesses from file.
  real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1)
  integer :: inconsistent = 0
  real :: dilate     ! The amount by which each layer is dilated to agree
                     ! with the bottom depth and free surface height, nondim.
  logical :: correct_thickness
  character(len=40)  :: mod = "initialize_thickness_from_file" ! This subroutine's name.
  character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path
  integer :: i, j, k, is, ie, js, je, nz

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, initialize_thickness_from_file: reading thickness", 5)

  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mod, "THICKNESS_FILE", thickness_file, &
                 "The name of the thickness file.", fail_if_missing=.true.)

  filename = trim(inputdir)//trim(thickness_file)
  call log_param(param_file, mod, "INPUTDIR/THICKNESS_FILE", filename)

  if (.not.file_exists(filename, G%Domain)) call GOLD_error(FATAL, &
         " initialize_thickness_from_file: Unable to open "//trim(filename))

  if (file_has_thickness) then
    call read_data(filename,"h",h(:,:,:),domain=G%Domain%mpp_domain)
  else
    call get_param(param_file, mod, "ADJUST_THICKNESS", correct_thickness, &
                 "If true, all mass below the bottom removed if the \n"//&
                 "topography is shallower than the thickness input file \n"//&
                 "would indicate.", default=.false.)

    call read_data(filename,"eta",eta(:,:,:),domain=G%Domain%mpp_domain)

    if (correct_thickness) then
      ! All mass below the bottom removed if the topography is shallower than
      ! the input file would indicate.  G%D is positive downward,
      ! eta is negative downward.
      do j=js,je ; do i=is,ie
        if (-eta(i,j,nz+1) > G%D(i,j) + 0.1) eta(i,j,nz+1) = -G%D(i,j)
      enddo ; enddo
    endif

    do k=nz,1,-1 ; do j=js,je ; do i=is,ie
      if (eta(i,j,K) < (eta(i,j,K+1) + G%Angstrom_z)) then
        eta(i,j,K) = eta(i,j,K+1) + G%Angstrom_z
        h(i,j,k) = G%Angstrom_z
      else
        h(i,j,k) = eta(i,j,K) - eta(i,j,K+1)
      endif
    enddo ; enddo ; enddo

  !  Check for consistency between the interface heights and topography.!
    if (correct_thickness) then
      do j=js,je ; do i=is,ie
        !   The whole column is dilated to accomodate deeper topography than
        ! the input file would indicate.
        if (-eta(i,j,nz+1) < G%D(i,j) - 0.1) then
          dilate = (eta(i,j,1)+G%D(i,j)) / (eta(i,j,1)-eta(i,j,nz+1))
          do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
        endif
      enddo ; enddo
    else
      do j=js,je ; do i=is,ie
        if (abs(eta(i,j,nz+1) + G%D(i,j)) > 1.0) inconsistent = inconsistent + 1
      enddo ; enddo
      call sum_across_PEs(inconsistent)

      if ((inconsistent > 0) .and. (is_root_pe())) then
        write(mesg,'("Thickness initial conditions are inconsistent ",'// &
                 '"with topography in ",I5," places.")') inconsistent
        call GOLD_error(WARNING, mesg)
      endif
    endif

  endif
end subroutine initialize_thickness_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_thickness_uniform(h, G, param_file)
  real, intent(out), dimension(NXMEM_,NYMEM_, NZ_) :: h
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file

! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!  This subroutine initializes the layer thicknesses to be uniform.
  real :: e0(SZK_(G)+1)   ! The resting interface heights, in m, usually !
                          ! negative because it is positive upward.      !
  real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface !
                          ! positive upward, in m.                       !
  real :: max_depth ! The maximum depths in m.
  integer :: i, j, k, is, ie, js, je, nz

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, initialize_thickness_uniform: setting thickness", 5)

  call get_param(param_file, "initialize_thickness_uniform", "MAXIMUM_DEPTH", &
                 max_depth, "The maximum depth of the ocean.", units="m", &
                 fail_if_missing=.true.)
  do k=1,nz
    e0(K) = -max_depth * real(k-1) / real(nz)
  enddo

  do j=js,je ; do i=is,ie
!    This sets the initial thickness (in m) of the layers.  The      !
!  thicknesses are set to insure that: 1.  each layer is at least an !
!  Angstrom thick, and 2.  the interfaces are where they should be   !
!  based on the resting depths and interface height perturbations,   !
!  as long at this doesn't interfere with 1.                         !
    eta1D(nz+1) = -1.0*G%D(i,j)
    do k=nz,1,-1
      eta1D(K) = e0(K)
      if (eta1D(K) < (eta1D(K+1) + G%Angstrom_z)) then
        eta1D(K) = eta1D(K+1) + G%Angstrom_z
        h(i,j,k) = G%Angstrom_z
      else
        h(i,j,k) = eta1D(K) - eta1D(K+1)
      endif
    enddo
  enddo ; enddo

end subroutine initialize_thickness_uniform
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_thickness_search
! search density space for location of layers
  call GOLD_error(FATAL,"  GOLD_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
end subroutine initialize_thickness_search
! -----------------------------------------------------------------------------

subroutine convert_thickness(h, G, param_file, tv)
  real, intent(inout), dimension(NXMEM_,NYMEM_, NZ_) :: h
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
  type(thermo_var_ptrs), intent(in) :: tv
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    p_top, p_bot
  real :: dz_geo(SZI_(G),SZJ_(G))      ! The change in geopotential height
                                       ! across a layer, in m2 s-2.
  real :: rho(SZI_(G))
  real :: I_gEarth
  logical :: Boussinesq
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: itt, max_itt

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  max_itt = 10
  Boussinesq = G%Boussinesq
  I_gEarth = 1.0 / G%g_Earth

  if (Boussinesq) then
    call GOLD_error(FATAL,"Not yet converting thickness with Boussinesq approx.")
  else
    if (associated(tv%eqn_of_state)) then
      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
        p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
      enddo ; enddo
      do k=1,nz
        do j=js,je
          do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, &
                                 is, ie-is+1, tv%eqn_of_state)
          do i=is,ie
            p_bot(i,j) = p_top(i,j) + G%g_Earth * h(i,j,k) * rho(i)
          enddo
        enddo

        do itt=1,max_itt
          call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, &
                                   0.0, G, tv%eqn_of_state, dz_geo)
          if (itt < max_itt) then ; do j=js,je
            call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, &
                                   is, ie-is+1, tv%eqn_of_state)
            ! Use Newton's method to correct the bottom value.
            !   The hydrostatic equation is linear to such a
            ! high degree that no bounds-checking is needed.
            do i=is,ie
              p_bot(i,j) = p_bot(i,j) + rho(i) * (G%g_Earth*h(i,j,k) - dz_geo(i,j))
            enddo
          enddo ; endif
        enddo

        do j=js,je ; do i=is,ie
          h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * G%kg_m2_to_H * I_gEarth
        enddo ; enddo
      enddo
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        h(i,j,k) = h(i,j,k) * G%Rlay(k) * G%kg_m2_to_H
      enddo ; enddo ; enddo
    endif
  endif

end subroutine convert_thickness

subroutine depress_surface(h, G, param_file, tv)
  real, dimension(NXMEM_,NYMEM_, NZ_), intent(inout) :: h
  type(ocean_grid_type),               intent(in)    :: G
  type(param_file_type),               intent(in)    :: param_file
  type(thermo_var_ptrs),               intent(in)    :: tv
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    eta_sfc  ! The free surface height that the model should use, in m.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    eta  ! The free surface height that the model should use, in m.
  real :: dilate  ! A ratio by which layers are dilated, nondim.
  real :: scale_factor ! A scaling factor for the eta_sfc values that are read
                       ! in, which can be used to change units, for example.
  character(len=40)  :: mod = "depress_surface" ! This subroutine's name.
  character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
  character(len=200) :: filename, eta_srf_var  ! Strings for file/path
  integer :: i, j, k, is, ie, js, je, nz
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

  ! Read the surface height (or pressure) from a file.

  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mod, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
                 "The initial condition file for the surface height.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
                 "The initial condition variable for the surface height.",&
                 default="SSH")
  filename = trim(inputdir)//trim(eta_srf_file)
  call log_param(param_file,  mod, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)

  call read_data(filename,eta_srf_var,eta_sfc,domain=G%Domain%mpp_domain)
  call get_param(param_file, mod, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
                 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into \n"//&
                 "units of m", units="variable", default=1.0)

  if (scale_factor /= 1.0) then ; do j=js,je ; do i=is,ie
    eta_sfc(i,j) = eta_sfc(i,j) * scale_factor
  enddo ; enddo ; endif

  ! Convert thicknesses to interface heights.
  call find_eta(h, tv, G%g_Earth, G, eta)
  
  do j=js,je ; do i=is,ie ; if (G%hmask(i,j) > 0.0) then
!    if (eta_sfc(i,j) < eta(i,j,nz+1)) then
      ! Issue a warning?
!    endif
    if (eta_sfc(i,j) > eta(i,j,1)) then
      ! Dilate the water column to agree, but only up to 10-fold.
      if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
        dilate = 10.0
        call GOLD_error(WARNING, "Free surface height dilation attempted "//&
               "to exceed 10-fold.", all_print=.true.)
      else
        dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
      endif
      do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
    elseif (eta(i,j,1) > eta_sfc(i,j)) then
      ! Remove any mass that is above the target free surface.
      do k=1,nz
        if (eta(i,j,K) <= eta_sfc(i,j)) exit
        if (eta(i,j,K+1) >= eta_sfc(i,j)) then
          h(i,j,k) = G%Angstrom
        else
          h(i,j,k) = max(G%Angstrom, h(i,j,k) * &
              (eta_sfc(i,j) - eta(i,j,K+1)) / (eta(i,j,K) - eta(i,j,K+1)) )
        endif
      enddo
    endif
  endif ; enddo ; enddo

end subroutine depress_surface

! -----------------------------------------------------------------------------
subroutine initialize_velocity_from_file(u, v, G, param_file)
  real, dimension(NXMEMQ_,NYMEM_, NZ_), intent(out) :: u
  real, dimension(NXMEM_,NYMEMQ_, NZ_), intent(out) :: v
  type(ocean_grid_type),                intent(in)  :: G
  type(param_file_type),                intent(in)  :: param_file
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine reads the initial velocity components from file
  character(len=40)  :: mod = "initialize_velocity_from_file" ! This subroutine's name.
  character(len=200) :: filename,velocity_file,inputdir ! Strings for file/path

  call GOLD_mesg("  GOLD_initialization.F90, initialize_velocity_from_file: reading u and v", 5)

  call get_param(param_file, mod, "VELOCITY_FILE", velocity_file, &
                 "The name of the velocity initial condition file.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)

  filename = trim(inputdir)//trim(velocity_file)
  call log_param(param_file, mod, "INPUTDIR/VELOCITY_FILE", filename)

  if (.not.file_exists(filename, G%Domain)) call GOLD_error(FATAL, &
         " initialize_velocity_from_file: Unable to open "//trim(filename))

  !  Read the velocities from a netcdf file.
  call read_data(filename,"u",u(:,:,:),domain=G%Domain%mpp_domain,position=EAST_FACE)
  call read_data(filename,"v",v(:,:,:),domain=G%Domain%mpp_domain,position=NORTH_FACE)

end subroutine initialize_velocity_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_velocity_zero(u, v, G, param_file)
  real, dimension(NXMEMQ_,NYMEM_, NZ_), intent(out) :: u
  real, dimension(NXMEM_,NYMEMQ_, NZ_), intent(out) :: v
  type(ocean_grid_type),                intent(in)  :: G
  type(param_file_type),                intent(in)  :: param_file
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to zero
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq


  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
    u(I,j,k) = 0.0
  enddo ; enddo ; enddo
  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
    v(i,J,k) = 0.0
  enddo ; enddo ; enddo

end subroutine initialize_velocity_zero
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_velocity_uniform(u, v, G, param_file)
  real, dimension(NXMEMQ_,NYMEM_, NZ_), intent(out) :: u
  real, dimension(NXMEM_,NYMEMQ_, NZ_), intent(out) :: v
  type(ocean_grid_type),                intent(in)  :: G
  type(param_file_type),                intent(in)  :: param_file
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to uniform
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  real    :: initial_u_const, initial_v_const
  character(len=200) :: mod = "initialize_velocity_uniform" ! This subroutine's name.
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq

  call get_param(param_file, mod, "INITIAL_U_CONST", initial_u_const, &
                 "A initial uniform value for the zonal flow.", &
                 units="m s-1", fail_if_missing=.true.)
  call get_param(param_file, mod, "INITIAL_V_CONST", initial_v_const, &
                 "A initial uniform value for the meridional flow.", &
                 units="m s-1", fail_if_missing=.true.)

  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
    u(I,j,k) = initial_u_const
  enddo ; enddo ; enddo
  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
    v(i,J,k) = initial_v_const
  enddo ; enddo ; enddo

end subroutine initialize_velocity_uniform
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_velocity_circular(u, v, G, param_file)
  real, dimension(NXMEMQ_,NYMEM_, NZ_), intent(out) :: u
  real, dimension(NXMEM_,NYMEMQ_, NZ_), intent(out) :: v
  type(ocean_grid_type),                intent(in)  :: G
  type(param_file_type),                intent(in)  :: param_file
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to be circular with
! no flow at edges of domain and center.
  character(len=200) :: mod = "initialize_velocity_circular"
  real :: south_lat, len_lat, west_lon, len_lon, circular_max_u
  real :: dpi, psi1, psi2
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq

  call get_param(param_file, mod, "SOUTHLAT", south_lat, &
                 "The southern latitude of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "LENLAT", len_lat, &
                 "The latitudinal length of the domain.", units="degrees", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "WESTLON", west_lon, &
                 "The western longitude of the domain.", units="degrees", &
                 default=0.0)
  call get_param(param_file, mod, "LENLON", len_lon, &
                 "The longitudinal length of the domain.", units="degrees", &
                 fail_if_missing=.true.)

  call get_param(param_file, mod, "CIRCULAR_MAX_U", circular_max_u, &
                 "The amplitude of zonal flow from which to scale the\n"// &
                 "circular stream function (m/s).", &
                 units="m s-1", default=0.)
  dpi=acos(0.0)*2.0 ! pi

  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
    psi1 = my_psi(I,j)
    psi2 = my_psi(I,j-1)
    u(I,j,k) = (psi1-psi2)/G%dy_u(I,j)! *(circular_max_u*len_lon/(2.0*dpi))
  enddo ; enddo ; enddo
  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
    psi1 = my_psi(i,J)
    psi2 = my_psi(i-1,J)
    v(i,J,k) = (psi2-psi1)/G%dx_v(i,J)! *(circular_max_u*len_lon/(2.0*dpi))
  enddo ; enddo ; enddo

  contains

  real function my_psi(ig,jg) ! in-line function
    integer :: ig, jg
    real :: x, y, r
    x = 2.0*(G%geolonq(ig,jg)-west_lon)/len_lon-1.0  ! -1<x<1
    y = 2.0*(G%geolatq(ig,jg)-south_lat)/len_lat-1.0 ! -1<y<1
    r = sqrt( x**2 + y**2 ) ! Circulat stream fn nis fn of radius only
    r = min(1.0,r) ! Flatten stream function in corners of box
    my_psi = 0.5*(1.0 - cos(dpi*r))
    my_psi = my_psi * (circular_max_u*len_lon*1e3/dpi) ! len_lon is in km
  end function my_psi

end subroutine initialize_velocity_circular
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_from_file(T, S, G, param_file)
  real, dimension(NXMEM_,NYMEM_, NZ_), intent(out) :: T, S
  type(ocean_grid_type),               intent(in)  :: G
  type(param_file_type),               intent(in)  :: param_file
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      from_file - .true. if the variables that are set here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  character(len=200) :: filename, ts_file, salt_file, inputdir ! Strings for file/path
  character(len=40)  :: mod = "initialize_temp_salt_from_file"
  character(len=64)  :: temp_var, salt_var

  call GOLD_mesg("  GOLD_initialization.F90, initialize_temp_salt_from_file: reading T and S", 5)

  call get_param(param_file, mod, "TS_FILE", ts_file, &
                 "The initial condition file for temperature.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)

  filename = trim(inputdir)//trim(ts_file)
  call log_param(param_file, mod, "INPUTDIR/TS_FILE", filename)
  if (.not.file_exists(filename, G%Domain)) call GOLD_error(FATAL, &
     " initialize_temp_salt_from_file: Unable to open "//trim(filename))

  call get_param(param_file, mod, "TEMP_IC_VAR", temp_var, &
                 "The initial condition variable for potential temperature.", &
                 default="PTEMP")
  call get_param(param_file, mod, "SALT_IC_VAR", salt_var, &
                 "The initial condition variable for salinity.", default="SALT")

! Read the temperatures and salinities from a netcdf file.           !
  call read_data(filename, temp_var, T(:,:,:), domain=G%Domain%mpp_domain)

  call get_param(param_file, mod, "SALT_FILE", salt_file, &
                 "The initial condition file for salinity.", default=trim(ts_file))
  filename = trim(inputdir)//trim(ts_file)
  if (.not.file_exists(filename, G%Domain)) call GOLD_error(FATAL, &
     " initialize_temp_salt_from_file: Unable to open "//trim(filename))

  call read_data(filename, salt_var, S(:,:,:), domain=G%Domain%mpp_domain)

end subroutine initialize_temp_salt_from_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_from_profile(T, S, G, param_file)
  real, dimension(NXMEM_,NYMEM_, NZ_), intent(out) :: T, S
  type(ocean_grid_type),               intent(in)  :: G
  type(param_file_type),               intent(in)  :: param_file
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      from_file - .true. if the variables that are set here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  real, dimension(SZK_(G)) :: T0, S0
  integer :: i, j, k
  character(len=200) :: filename, ts_file, inputdir ! Strings for file/path
  character(len=40)  :: mod = "initialize_temp_salt_from_profile"

  call GOLD_mesg("  GOLD_initialization.F90, initialize_temp_salt_from_file: reading T and S", 5)

  call get_param(param_file, mod, "TS_FILE", ts_file, &
                 "The file with the reference profiles for temperature \n"//&
                 "and salinity.", fail_if_missing=.true.)
  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(ts_file)
  call log_param(param_file, mod, "INPUTDIR/TS_FILE", filename)
  if (.not.file_exists(filename)) call GOLD_error(FATAL, &
     " initialize_temp_salt_from_profile: Unable to open "//trim(filename))

! Read the temperatures and salinities from a netcdf file.           !
  call read_data(filename,"PTEMP",T0(:),domain=G%Domain%mpp_domain)
  call read_data(filename,"SALT", S0(:),domain=G%Domain%mpp_domain)

  do k=1,G%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
    T(i,j,k) = T0(k) ; S(i,j,k) = S0(k)
  enddo ; enddo ; enddo

end subroutine initialize_temp_salt_from_profile
! -----------------------------------------------------------------------------


! -----------------------------------------------------------------------------
subroutine initialize_temp_salt_fit(T, S, G, param_file, eqn_of_state, P_Ref)
  real, dimension(NXMEM_,NYMEM_, NZ_), intent(out) :: T, S
  type(ocean_grid_type),               intent(in)  :: G
  type(param_file_type),               intent(in)  :: param_file
  type(EOS_type),                      pointer     :: eqn_of_state
  real,                                intent(in)  :: P_Ref
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      eqn_of_state - integer that selects the equatio of state
!  (in)      P_Ref - The coordinate-density reference pressure in Pa.
  real :: T0(SZK_(G)), S0(SZK_(G))
  real :: T_Ref         ! Reference Temperature
  real :: S_Ref         ! Reference Salinity
  real :: pres(SZK_(G))      ! An array of the reference pressure in Pa.
  real :: drho_dT(SZK_(G))   ! Derivative of density with temperature in kg m-3 K-1.                              !
  real :: drho_dS(SZK_(G))   ! Derivative of density with salinity in kg m-3 PSU-1.                             !
  real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 in kg m-3.
  character(len=40)  :: mod = "initialize_temp_salt_fit" ! This subroutine's name.
  integer :: i, j, k, itt, nz
  nz = G%ke

  call GOLD_mesg("  GOLD_initialization.F90, initialize_temp_salt_fit: "//&
                 "  setting T and S", 5)

  call get_param(param_file, mod, "T_REF", T_Ref, &
                 "A reference temperature used in initialization.", &
                 units="degC", fail_if_missing=.true.)
  call get_param(param_file, mod, "S_REF", S_Ref, &
                 "A reference salinity used in initialization.", units="PSU", &
                 default=35.0)
  do k=1,nz
    pres(k) = P_Ref ; S0(k) = S_Ref
  enddo
  T0(1) = T_Ref

  call calculate_density(T0(1),S0(1),pres(1),rho_guess(1),1,1,eqn_of_state)
  call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,eqn_of_state)

! A first guess of the layers' temperatures.                         !
  do k=nz,1,-1
    T0(k) = T0(1) + (G%Rlay(k) - rho_guess(1)) / drho_dT(1)
  enddo

! Refine the guesses for each layer.                                 !
  do itt=1,6
    call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state)
    call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state)
    do k=1,nz
      T0(k) = T0(k) + (G%Rlay(k) - rho_guess(k)) / drho_dT(k)
    enddo
  enddo

  do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
    T(i,j,k) = T0(k) ; S(i,j,k) = S0(k)
  enddo ; enddo ; enddo

end subroutine initialize_temp_salt_fit
! -----------------------------------------------------------------------------


! -----------------------------------------------------------------------------
subroutine initialize_sponges_file(G, use_temperature, tv, param_file, CSp)
  type(ocean_grid_type), intent(in) :: G
  logical, intent(in) :: use_temperature
  type(thermo_var_ptrs), intent(in) :: tv
  type(param_file_type), intent(in) :: param_file
  type(sponge_CS),       pointer    :: CSp
!   This subroutine sets the inverse restoration time (Idamp), and   !
! the values towards which the interface heights and an arbitrary    !
! number of tracers should be restored within each sponge. The       !
! interface height is always subject to damping, and must always be  !
! the first registered field.                                        !

! Arguments: from_file - .true. if the variables that are used here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read for all fields
!                       except the inverse damping rate.
!  (in)      damp_file - The name of the file from which to read the
!                        inverse damping rate.
!  (in)      G - The ocean's grid structure.
!  (in)      use_temperature - If true, T & S are state variables.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in/out)  CSp - A pointer that is set to point to the control structure
!                  for this module

  real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! The target interface heights, in m.
  real, dimension (SZI_(G),SZJ_(G),SZK_(G)) :: &
    tmp, tmp2 ! A temporary array for tracers.
  real, dimension (SZI_(G),SZJ_(G)) :: &
    tmp_2d ! A temporary array for tracers.

  real :: Idamp(SZI_(G),SZJ_(G))    ! The inverse damping rate, in s-1.
  real :: pres(SZI_(G))     ! An array of the reference pressure, in Pa.

  integer :: i, j, k, is, ie, js, je, nz
  character(len=40) :: potemp_var, salin_var, Idamp_var, eta_var
  character(len=40) :: mod = "initialize_sponges_file"
  character(len=200) :: damping_file, state_file  ! Strings for filenames
  character(len=200) :: filename, inputdir ! Strings for file/path and path.
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

  pres(:) = 0.0 ; eta(:,:,:) = 0.0 ; tmp(:,:,:) = 0.0 ; Idamp(:,:) = 0.0

  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mod, "SPONGE_DAMPING_FILE", damping_file, &
                 "The name of the file with the sponge damping rates.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "SPONGE_STATE_FILE", state_file, &
                 "The name of the file with the state to damp toward.", &
                 default=damping_file)

  call get_param(param_file, mod, "SPONGE_PTEMP_VAR", potemp_var, &
                 "The name of the potential temperature variable in \n"//&
                 "SPONGE_STATE_FILE.", default="PTEMP")
  call get_param(param_file, mod, "SPONGE_SALT_VAR", salin_var, &
                 "The name of the salinity variable in \n"//&
                 "SPONGE_STATE_FILE.", default="SALT")
  call get_param(param_file, mod, "SPONGE_ETA_VAR", eta_var, &
                 "The name of the interface height variable in \n"//&
                 "SPONGE_STATE_FILE.", default="ETA")
  call get_param(param_file, mod, "SPONGE_IDAMP_VAR", Idamp_var, &
                 "The name of the inverse damping rate variable in \n"//&
                 "SPONGE_DAMPING_FILE.", default="IDAMP")

  filename = trim(inputdir)//trim(damping_file)
  call log_param(param_file, mod, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
  if (.not.file_exists(filename, G%Domain)) &
    call GOLD_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))


  call read_data(filename,"Idamp",Idamp(:,:), domain=G%Domain%mpp_domain)

! Now register all of the fields which are damped in the sponge.     !
! By default, momentum is advected vertically within the sponge, but !
! momentum is typically not damped within the sponge.                !

  filename = trim(inputdir)//trim(state_file)
  call log_param(param_file, mod, "INPUTDIR/SPONGE_STATE_FILE", filename)
  if (.not.file_exists(filename, G%Domain)) &
    call GOLD_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))


!  The first call to set_up_sponge_field is for the interface height.!
  call read_data(filename, eta_var, eta(:,:,:), domain=G%Domain%mpp_domain)

  do j=js,je ; do i=is,ie
    eta(i,j,nz+1) = -G%D(i,j)
  enddo ; enddo
  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
    if (eta(i,j,K) < (eta(i,j,K+1) + G%Angstrom_z)) &
      eta(i,j,K) = eta(i,j,K+1) + G%Angstrom_z
  enddo ; enddo ; enddo
! Set the inverse damping rates so that the model will know where to !
! apply the sponges, along with the interface heights.               !
  call initialize_sponge(Idamp, eta, G, param_file, CSp)

!   Now register all of the tracer fields which are damped in the    !
! sponge. By default, momentum is advected vertically within the     !
! sponge, but momentum is typically not damped within the sponge.    !

  if ( G%nkml>0 ) then
!   This call to set_up_sponge_ML_density registers the target values of the
! mixed layer density, which is used in determining which layers can be
! inflated without causing static instabilities.
    do i=is-1,ie ; pres(i) = tv%P_Ref ; enddo

    call read_data(filename, potemp_var, tmp(:,:,:), domain=G%Domain%mpp_domain)
    call read_data(filename, salin_var, tmp2(:,:,:), domain=G%Domain%mpp_domain)

    do j=js,je
      call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), &
                             is, ie-is+1, tv%eqn_of_state)
    enddo

    call set_up_sponge_ML_density(tmp_2d, CSp)
  endif

!  The remaining calls to set_up_sponge_field can be in any order.   !
  if ( use_temperature ) then
    call read_data(filename, potemp_var, tmp(:,:,:), domain=G%Domain%mpp_domain)
    call set_up_sponge_field(tmp, tv%T, nz, CSp)
    call read_data(filename, salin_var, tmp(:,:,:), domain=G%Domain%mpp_domain)
    call set_up_sponge_field(tmp, tv%S, nz, CSp)
  endif


end subroutine initialize_sponges_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_Open_Bdry_Conds(OBC, tv, G, param_file, advect_tracer_CSp)
  type(ocean_OBC_type),  pointer    :: OBC
  type(thermo_var_ptrs), intent(in) :: tv
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
  type(advect_tracer_CS), pointer   :: advect_tracer_CSp
!   This subroutine sets the properties of flow at open boundary conditions.
! This particular example is for the DOME inflow describe in Legg et al. 2006.

! Arguments: OBC - This open boundary condition type specifies whether, where,
!                  and what open boundary conditions are used.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

  logical :: any_OBC        ! Set to true if any points in this subdomain use
                            ! open boundary conditions.
  logical, pointer, dimension(:,:) :: &
    OBC_mask_u => NULL(), & ! These arrays are true at zonal or meridional
    OBC_mask_v => NULL()    ! velocity points that have prescribed open boundary
                            ! conditions.
  real, pointer, dimension(:,:,:) :: &
    OBC_T_u => NULL(), &    ! These arrays should be allocated and set to
    OBC_T_v => NULL(), &    ! specify the values of T and S that should come
    OBC_S_u => NULL(), &    ! in through u- and v- points through the open
    OBC_S_v => NULL()       ! boundary conditions, in C and psu.
  logical :: apply_OBC_u = .false., apply_OBC_v = .false.
  ! The following variables are used to set the target temperature and salinity.
  real :: T0(SZK_(G)), S0(SZK_(G))
  real :: pres(SZK_(G))      ! An array of the reference pressure in Pa.
  real :: drho_dT(SZK_(G))   ! Derivative of density with temperature in kg m-3 K-1.                              !
  real :: drho_dS(SZK_(G))   ! Derivative of density with salinity in kg m-3 PSU-1.                             !
  real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 in kg m-3.
  character(len=40) :: mod = "set_Open_Bdry_Conds"
  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, yhalo
  integer :: Isdq, Iedq, Jsdq, Jedq

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq
  yhalo = G%jsc-G%jsd

  call get_param(param_file, mod, "APPLY_OBC_U", apply_OBC_u, default=.false.)
  call get_param(param_file, mod, "APPLY_OBC_V", apply_OBC_v, default=.false.)

  if (apply_OBC_u) then
    ! Determine where u points are applied.
    allocate(OBC_mask_u(Isdq:Iedq,jsd:jed)) ; OBC_mask_u(:,:) = .false.
    any_OBC = .false.
    do j=jsd,jed ; do I=Isdq,Iedq
    ! if (SOME_TEST_FOR_U_OPEN_BCS) then
    !   OBC_mask_u(I,j) = .true. ; any_OBC = .true.
    ! endif
    enddo ; enddo
    if (.not.any_OBC) then
      ! This processor has no u points at which open boundary conditions are
      ! to be applied.
      apply_OBC_u = .false.
      deallocate(OBC_mask_u)
    endif
  endif
  if (apply_OBC_v) then
    ! Determine where v points are applied.
    allocate(OBC_mask_v(isd:ied,Jsdq:Jedq)) ; OBC_mask_v(:,:) = .false.
    any_OBC = .false.
    do J=Jsdq,Jedq ; do i=isd,ied
    ! if (SOME_TEST_FOR_V_OPEN_BCS) then
    !   OBC_mask_v(i,J) = .true. ; any_OBC = .true.
    ! endif
    enddo ; enddo
    if (.not.any_OBC) then
      ! This processor has no v points at which open boundary conditions are
      ! to be applied.
      apply_OBC_v = .false.
      deallocate(OBC_mask_v)
    endif
  endif

  if (.not.(apply_OBC_u .or. apply_OBC_v)) return

  if (.not.associated(OBC)) allocate(OBC)

  if (apply_OBC_u) then
    OBC%apply_OBC_u = .true.
    OBC%OBC_mask_u => OBC_mask_u
    allocate(OBC%u(Isdq:Iedq,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0
    allocate(OBC%uh(Isdq:Iedq,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0
    allocate(OBC%OBC_kind_u(Isdq:Iedq,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE
    do j=jsd,jed ; do I=Isdq,Iedq
      if (OBC%OBC_mask_u(I,j)) OBC%OBC_kind_u(I,j) = OBC_SIMPLE
    enddo ; enddo
  endif
  if (apply_OBC_v) then
    OBC%apply_OBC_v = .true.
    OBC%OBC_mask_v => OBC_mask_v
    allocate(OBC%v(isd:ied,Jsdq:Jedq,nz)) ; OBC%v(:,:,:) = 0.0
    allocate(OBC%vh(isd:ied,Jsdq:Jedq,nz)) ; OBC%vh(:,:,:) = 0.0
    allocate(OBC%OBC_kind_v(isd:ied,Jsdq:Jedq)) ; OBC%OBC_kind_v(:,:) = OBC_NONE
    do J=Jsdq,Jedq ; do i=isd,ied
      if (OBC%OBC_mask_v(i,J)) OBC%OBC_kind_v(i,J) = OBC_SIMPLE
    enddo ; enddo
  endif

  if (apply_OBC_v) then
    do k=1,nz ; do J=Jsd,Jed ; do i=isd,ied
      if (OBC_mask_v(i,J)) then
        ! An appropriate expression for the meridional inflow velocities and
        ! transports should go here.
        OBC%vh(i,J,k) = 0.0 * G%m_to_H ; OBC%v(i,J,k) = 0.0
      else
        OBC%vh(i,J,k) = 0.0 ; OBC%v(i,J,k) = 0.0
      endif
    enddo ; enddo ; enddo
  endif

  if (apply_OBC_u) then
    do k=1,nz ; do j=jsd,jed ; do I=Isdq,Iedq
      if (OBC_mask_u(I,j)) then
        ! An appropriate expression for the zonal inflow velocities and
        ! transports should go here.
        OBC%uh(I,j,k) = 0.0 * G%m_to_H ; OBC%u(I,j,k) = 0.0
      else
        OBC%uh(I,j,k) = 0.0 ; OBC%u(I,j,k) = 0.0
      endif
    enddo ; enddo ; enddo
  endif

  !   The inflow values of temperature and salinity also need to be set here if
  ! these variables are used.  The following code is just a naive example.
  if (apply_OBC_u .or. apply_OBC_v) then
    if (associated(tv%S)) then
      ! In this example, all S inflows have values of 35 psu.
      call add_tracer_OBC_values("S", advect_tracer_CSp, OBC_inflow=35.0)
    endif
    if (associated(tv%T)) then
      ! In this example, the T values are set to be consistent with the layer
      ! target density and a salinity of 35 psu.  This code is taken from
      !  initialize_temp_sal.
      pres(:) = tv%P_Ref ; S0(:) = 35.0 ; T0(1) = 25.0
      call calculate_density(T0(1),S0(1),pres(1),rho_guess(1),1,1,tv%eqn_of_state)
      call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,tv%eqn_of_state)

      do k=1,nz ; T0(k) = T0(1) + (G%Rlay(k)-rho_guess(1)) / drho_dT(1) ; enddo
      do itt=1,6
        call calculate_density(T0,S0,pres,rho_guess,1,nz,tv%eqn_of_state)
        call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,tv%eqn_of_state)
        do k=1,nz ; T0(k) = T0(k) + (G%Rlay(k)-rho_guess(k)) / drho_dT(k) ; enddo
      enddo

      if (apply_OBC_u) then
        allocate(OBC_T_u(Isdq:Iedq,jsd:jed,nz))
        do k=1,nz ; do j=jsd,jed ; do I=Isdq,Iedq
          OBC_T_u(I,j,k) = T0(k)
        enddo ; enddo ; enddo
      endif
      if (apply_OBC_v) then
        allocate(OBC_T_v(isd:ied,Jsdq:Jedq,nz))
        do k=1,nz ; do J=Jsdq,Jedq ; do i=isd,ied
          OBC_T_v(i,J,k) = T0(k)
        enddo ; enddo ; enddo
      endif
      call add_tracer_OBC_values("T", advect_tracer_CSp, OBC_in_u=OBC_T_u, &
                                                         OBC_in_v=OBC_T_v)
    endif
  endif

end subroutine set_Open_Bdry_Conds
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, advect_tracer_CSp)
  type(ocean_OBC_type),                pointer    :: OBC
  type(thermo_var_ptrs),               intent(inout) :: tv
  real, dimension(NXMEM_,NYMEM_, NZ_), intent(inout) :: h
  type(ocean_grid_type),               intent(inout) :: G
  type(param_file_type),               intent(in) :: PF
  type(advect_tracer_CS),              pointer    :: advect_tracer_CSp
!   This subroutine sets the initial definitions of the characteristic open boundary
!   conditions. Written by Mehmet Ilicak

! Arguments: OBC - This open boundary condition type specifies whether, where,
!                  and what open boundary conditions are used.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (in)      G - The ocean's grid structure.
!  (in)      PF - A structure indicating the open file to parse for
!                         model parameter values.

  logical :: any_OBC        ! Set to true if any points in this subdomain use
                            ! open boundary conditions.

  logical :: apply_OBC_u_flather_east = .false., apply_OBC_u_flather_west = .false.
  logical :: apply_OBC_v_flather_north = .false., apply_OBC_v_flather_south = .false.  
  logical :: read_OBC_eta = .false.
  logical :: read_OBC_uv = .false.
  logical :: read_OBC_TS = .false.

  integer :: isd_global, jsd_global
  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, yhalo
  integer :: isd_off, jsd_off
  integer :: Isdq, Iedq, Jsdq, Jedq
  integer :: east_boundary, west_boundary, north_boundary, south_boundary
  character(len=40)  :: mod = "set_Flather_Bdry_Conds" ! This subroutine's name.
  character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path

  real :: temp_u(G%domain%nxtot+1,G%domain%nytot)
  real :: temp_v(G%domain%nxtot,G%domain%nytot+1)

  real, pointer, dimension(:,:,:) :: &
    OBC_T_u => NULL(), &    ! These arrays should be allocated and set to
    OBC_T_v => NULL(), &    ! specify the values of T and S that should come
    OBC_S_u => NULL(), & 
    OBC_S_v => NULL()     

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq
  yhalo = G%jsc-G%jsd 
  
  isd_global = G%isd_global
  jsd_global = G%jsd_global
  call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,&
                 "Apply a Flather open boundary condition on the eastern \n"//&
                 "side of the global domain", default=.false.)
  call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,&
                 "Apply a Flather open boundary condition on the western \n"//&
                 "side of the global domain", default=.false.)
  call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,&
                 "Apply a Flather open boundary condition on the northern \n"//&
                 "side of the global domain", default=.false.)
  call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,&
                 "Apply a Flather open boundary condition on the southern \n"//&
                 "side of the global domain", default=.false.)

  if (.not.(apply_OBC_u_flather_east  .or. apply_OBC_u_flather_west .or. &            
            apply_OBC_v_flather_north .or. apply_OBC_v_flather_south)) return
  
  if (.not.associated(OBC)) allocate(OBC)
       
  OBC%apply_OBC_u_flather_east = apply_OBC_u_flather_east
  OBC%apply_OBC_u_flather_west = apply_OBC_u_flather_west 
  OBC%apply_OBC_v_flather_north = apply_OBC_v_flather_north 
  OBC%apply_OBC_v_flather_south = apply_OBC_v_flather_south

  call get_param(PF, mod, "READ_OBC_UV", read_OBC_uv, &
                 "If true, read the values for the velocity open boundary \n"//&
                 "conditions from the file specified by OBC_FILE.", &
                 default=.false.)
  call get_param(PF, mod, "READ_OBC_ETA", read_OBC_eta, &
                 "If true, read the values for the sea surface height \n"//&
                 "open boundary conditions from the file specified by \n"//&
                 "OBC_FILE.", default=.false.)
  call get_param(PF, mod, "READ_OBC_TS", read_OBC_TS, &
                 "If true, read the values for the temperature and \n"//&
                 "salinity open boundary conditions from the file \n"//&
                 "specified by OBC_FILE.", default=.false.)
  if (read_OBC_uv .or. read_OBC_eta .or. read_OBC_TS) then
    call get_param(PF, mod, "OBC_FILE", OBC_file, &
                 "The file from which the appropriate open boundary \n"//&
                 "condition values are read.", default="GOLD_OBC_FILE.nc")
    call get_param(PF, mod, "INPUTDIR", inputdir, default=".")
    inputdir = slasher(inputdir)
    filename = trim(inputdir)//trim(OBC_file)
    call log_param(PF, mod, "INPUTDIR/OBC_FILE", filename)
  endif

  if (G%symmetric) then
    east_boundary = G%domain%nxtot+G%domain%nx_halo
    west_boundary = G%domain%nx_halo
    north_boundary = G%domain%nytot+G%domain%ny_halo
    south_boundary = G%domain%ny_halo
  else
    east_boundary = G%domain%nxtot+G%domain%nx_halo-1
    west_boundary = G%domain%nx_halo+1
    north_boundary = G%domain%nytot+G%domain%ny_halo-1
    south_boundary = G%domain%ny_halo+1
  endif

  if (.not.associated(OBC%OBC_mask_u)) then
    allocate(OBC%OBC_mask_u(Isdq:Iedq,jsd:jed)) ; OBC%OBC_mask_u(:,:) = .false.
  endif
  if (.not.associated(OBC%OBC_kind_u)) then
    allocate(OBC%OBC_kind_u(Isdq:Iedq,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE
  endif
  if (.not.associated(OBC%OBC_mask_v)) then
    allocate(OBC%OBC_mask_v(isd:ied,Jsdq:Jedq)) ; OBC%OBC_mask_v(:,:) = .false.
  endif
  if (.not.associated(OBC%OBC_kind_v)) then
    allocate(OBC%OBC_kind_v(isd:ied,Jsdq:Jedq)) ; OBC%OBC_kind_v(:,:) = OBC_NONE
  endif

  if (.not.associated(OBC%vbt_outer)) then
    allocate(OBC%vbt_outer(isd:ied,Jsdq:Jedq)) ; OBC%vbt_outer(:,:) = 0.0
  endif

  if (.not.associated(OBC%ubt_outer)) then
    allocate(OBC%ubt_outer(Isdq:Iedq,jsd:jed)) ; OBC%ubt_outer(:,:) = 0.0
  endif

  if (.not.associated(OBC%eta_outer_u)) then
    allocate(OBC%eta_outer_u(Isdq:Iedq,jsd:jed)) ; OBC%eta_outer_u(:,:) = 0.0
  endif

  if (.not.associated(OBC%eta_outer_v)) then
    allocate(OBC%eta_outer_v(isd:ied,Jsdq:Jedq)) ; OBC%eta_outer_v(:,:) = 0.0
  endif
  
  if (read_OBC_uv) then
    call read_data(filename, 'ubt', OBC%ubt_outer, &
                   domain=G%Domain%mpp_domain, position=EAST_FACE)
    call read_data(filename, 'vbt', OBC%vbt_outer, &
                   domain=G%Domain%mpp_domain, position=NORTH_FACE)
  endif

  if (read_OBC_eta) then
    call read_data(filename, 'eta_outer_u', OBC%eta_outer_u, &
                   domain=G%Domain%mpp_domain, position=EAST_FACE)	
    call read_data(filename, 'eta_outer_v', OBC%eta_outer_v, &
                   domain=G%Domain%mpp_domain, position=NORTH_FACE)
  endif

  call pass_vector(OBC%eta_outer_u,OBC%eta_outer_v,G%Domain, To_All+SCALAR_PAIR, CGRID_NE)
  call pass_vector(OBC%ubt_outer,OBC%vbt_outer,G%Domain)

  ! This code should be modified to allow OBCs to be applied anywhere.

  if (apply_OBC_u_flather_east) then
    ! Determine where u points are applied at east side 
    do j=jsd,jed ; do I=Isdq,Iedq
      if ((I+isd_global-isd) .eq. east_boundary) then !eastern side
        OBC%OBC_mask_u(I,j) = .true.
        OBC%OBC_kind_u(I,j) = OBC_FLATHER_E
        if ((i+1>isd) .and. (i+1<ied) .and. (J>Jsdq) .and. (J<Jedq)) then
          OBC%OBC_mask_v(i+1,J) = .true.
          if (OBC%OBC_kind_v(i+1,J) == OBC_NONE) OBC%OBC_kind_v(i+1,J) = OBC_FLATHER_E
        endif
        if ((i+1>isd) .and. (i+1<ied) .and. (J-1>Jsdq) .and. (J-1<Jedq)) then
          OBC%OBC_mask_v(i+1,J-1) = .true.
          if (OBC%OBC_kind_v(i+1,J-1) == OBC_NONE) OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER_E
        endif
      endif
    enddo ; enddo
  endif

  if (apply_OBC_u_flather_west) then
    ! Determine where u points are applied at west side 
    do j=jsd,jed ; do I=Isdq,Iedq
      if ((I+isd_global-isd) .eq. west_boundary) then !western side
        OBC%OBC_mask_u(I,j) = .true.
        OBC%OBC_kind_u(I,j) = OBC_FLATHER_W
        if ((i>isd) .and. (i<ied) .and. (J>Jsdq) .and. (J<Jedq)) then
          OBC%OBC_mask_v(i,J) = .true.
          if (OBC%OBC_kind_v(i,J) == OBC_NONE) OBC%OBC_kind_v(i,J) = OBC_FLATHER_W
        endif
        if ((i>isd) .and. (i<ied) .and. (J-1>Jsdq) .and. (J-1<Jedq)) then
          OBC%OBC_mask_v(i,J-1) = .true.
          if (OBC%OBC_kind_v(i,J-1) == OBC_NONE) OBC%OBC_kind_v(i,J-1) = OBC_FLATHER_W
        endif
      endif
    enddo ; enddo
  endif


  if (apply_OBC_v_flather_north) then
    ! Determine where v points are applied at north side 
    do J=Jsdq,Jedq ; do i=isd,ied
      if ((J+jsd_global-jsd) .eq. north_boundary) then         !northern side
        OBC%OBC_mask_v(i,J) = .true.
        OBC%OBC_kind_v(i,J) = OBC_FLATHER_N
        if ((I>Isdq) .and. (I<Iedq) .and. (j+1>jsd) .and. (j+1<jed)) then
          OBC%OBC_mask_u(I,j+1) = .true.
          if (OBC%OBC_kind_u(I,j+1) == OBC_NONE) OBC%OBC_kind_u(I,j+1) = OBC_FLATHER_N
        endif
        if ((I-1>Isdq) .and. (I-1<Iedq) .and. (j+1>jsd) .and. (j+1<jed)) then
          OBC%OBC_mask_u(I-1,j+1) = .true.
          if (OBC%OBC_kind_u(I-1,j+1) == OBC_NONE) OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER_N
        endif
     endif
    enddo ; enddo
  endif
  
  if (apply_OBC_v_flather_south) then
    ! Determine where v points are applied at south side 
    do J=Jsdq,Jedq ; do i=isd,ied
      if ((J+jsd_global-jsd) .eq. south_boundary) then         !southern side
        OBC%OBC_mask_v(i,J) = .true.
        OBC%OBC_kind_v(i,J) = OBC_FLATHER_S
        if ((I>Isdq) .and. (I<Iedq) .and. (j>jsd) .and. (j<jed)) then
          OBC%OBC_mask_u(I,j) = .true.
          if (OBC%OBC_kind_u(I,j) == OBC_NONE) OBC%OBC_kind_u(I,j) = OBC_FLATHER_S
        endif
        if ((I-1>Isdq) .and. (I-1<Iedq) .and. (j>jsd) .and. (j<jed)) then
          OBC%OBC_mask_u(I-1,j) = .true.
          if (OBC%OBC_kind_u(I-1,j) == OBC_NONE) OBC%OBC_kind_u(I-1,j) = OBC_FLATHER_S
        endif
      endif
    enddo ; enddo
  endif

  !   If there are no OBC points on this PE, there is no reason to keep the OBC
  ! type, and it could be deallocated.


  ! Define radiation coefficients r[xy]_old_[uvh] as needed.  For now, there are
  ! no radiation conditions applied to the thicknesses, since the thicknesses
  ! might not be physically motivated.  Instead, sponges should be used to
  ! enforce the near-boundary layer structure.
  if (apply_OBC_u_flather_west .or. apply_OBC_u_flather_east) then
    allocate(OBC%rx_old_u(Isdq:Iedq,jsd:jed,nz)) ; OBC%rx_old_u(:,:,:) = 0.0
 !   allocate(OBC%rx_old_h(Isd:Ied,jsd:jed,nz))   ; OBC%rx_old_h(:,:,:) = 0.0
  endif
  if (apply_OBC_v_flather_south .or. apply_OBC_v_flather_north) then
    allocate(OBC%ry_old_v(isd:ied,Jsdq:Jedq,nz)) ; OBC%ry_old_v(:,:,:) = 0.0
 !   allocate(OBC%ry_old_h(isd:ied,Jsd:Jed,nz))   ; OBC%ry_old_h(:,:,:) = 0.0
  endif


  if (associated(tv%T)) then
    allocate(OBC_T_u(Isdq:Iedq,jsd:jed,nz)) ; OBC_T_u(:,:,:) = 0.0
    allocate(OBC_S_u(Isdq:Iedq,jsd:jed,nz)) ; OBC_S_u(:,:,:) = 0.0
    allocate(OBC_T_v(isd:ied,Jsdq:Jedq,nz)) ; OBC_T_v(:,:,:) = 0.0
    allocate(OBC_S_v(isd:ied,Jsdq:Jedq,nz)) ; OBC_S_v(:,:,:) = 0.0

    if (read_OBC_TS) then
      call read_data(filename, 'OBC_T_u', OBC_T_u, &
                     domain=G%Domain%mpp_domain, position=EAST_FACE)
      call read_data(filename, 'OBC_S_u', OBC_S_u, &
                     domain=G%Domain%mpp_domain, position=EAST_FACE)

      call read_data(filename, 'OBC_T_v', OBC_T_v, &
                     domain=G%Domain%mpp_domain, position=NORTH_FACE)
      call read_data(filename, 'OBC_S_v', OBC_S_v, &
                     domain=G%Domain%mpp_domain, position=NORTH_FACE)
    else
      call pass_var(tv%T, G%Domain)
      call pass_var(tv%S, G%Domain)
      do k=1,nz ; do j=js,je ; do I=is-1,ie
        if (OBC%OBC_mask_u(I,j)) then
          if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) then
            OBC_T_u(I,j,k) = tv%T(i,j,k)
            OBC_S_u(I,j,k) = tv%S(i,j,k)
          elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then
            OBC_T_u(I,j,k) = tv%T(i+1,j,k)
            OBC_S_u(I,j,k) = tv%S(i+1,j,k)
          elseif (G%hmask(i,j) + G%hmask(i+1,j) > 0) then
            OBC_T_u(I,j,k) = (G%hmask(i,j)*tv%T(i,j,k) + G%hmask(i+1,j)*tv%T(i+1,j,k)) / &
                             (G%hmask(i,j) + G%hmask(i+1,j))
            OBC_S_u(I,j,k) = (G%hmask(i,j)*tv%S(i,j,k) + G%hmask(i+1,j)*tv%S(i+1,j,k)) / &
                             (G%hmask(i,j) + G%hmask(i+1,j))
          else ! This probably shouldn't happen or maybe it doesn't matter?
            OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k))
            OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k))
          endif
        else
          OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k))
          OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k))
        endif
      enddo; enddo ; enddo

      do k=1,nz ; do J=js-1,je ; do i=is,ie
        if (OBC%OBC_mask_v(i,J)) then
          if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) then
            OBC_T_v(i,J,k) = tv%T(i,j,k)
            OBC_S_v(i,J,k) = tv%S(i,j,k)
          elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then
            OBC_T_v(i,J,k) = tv%T(i,j+1,k)
            OBC_S_v(i,J,k) = tv%S(i,j+1,k)
          elseif (G%hmask(i,j) + G%hmask(i,j+1) > 0) then
            OBC_T_v(i,J,k) = (G%hmask(i,j)*tv%T(i,j,k) + G%hmask(i,j+1)*tv%T(i,j+1,k)) / &
                             (G%hmask(i,j) + G%hmask(i,j+1))
            OBC_S_v(i,J,k) = (G%hmask(i,j)*tv%S(i,j,k) + G%hmask(i,j+1)*tv%S(i,j+1,k)) / &
                             (G%hmask(i,j) + G%hmask(i,j+1))
          else ! This probably shouldn't happen or maybe it doesn't matter?
            OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k))
            OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k))
          endif
        else
          OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k))
          OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k))
        endif
      enddo; enddo ; enddo
    endif

    call pass_vector(OBC_T_u, OBC_T_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)
    call pass_vector(OBC_S_u, OBC_S_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)

    call add_tracer_OBC_values("T", advect_tracer_CSp, OBC_in_u=OBC_T_u, &
                                                       OBC_in_v=OBC_T_v)
    call add_tracer_OBC_values("S", advect_tracer_CSp, OBC_in_u=OBC_S_u, &
                                                       OBC_in_v=OBC_S_v)
    do k=1,nz ; do j=js,je ; do I=is-1,ie
      if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) then
        tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k)
      elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then
        tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k)
      endif
    enddo ; enddo ; enddo
    do k=1,nz ; do J=js-1,je ; do i=is,ie
      if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) then
        tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k)
      elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then
        tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k)
      endif
    enddo ; enddo ; enddo
  endif

  do k=1,nz ; do j=js,je ; do I=is-1,ie
    if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) h(i+1,j,k) = h(i,j,k)
    if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) h(i,j,k) = h(i+1,j,k)
  enddo ; enddo ; enddo
  do k=1,nz ; do J=js-1,je ; do i=is,ie
    if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) h(i,j+1,k) = h(i,j,k)
    if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) h(i,j,k) = h(i,j+1,k)
  enddo ; enddo ; enddo

end subroutine set_Flather_Bdry_Conds   
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine reset_face_lengths_named(G, param_file, name)
  type(ocean_grid_type), intent(inout) :: G
  type(param_file_type), intent(in) :: param_file
  character(len=*),      intent(in) :: name
!   This subroutine sets the open face lengths at selected points to restrict
! passages to their observed widths.

! Arguments: G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      name - The name for the set of face lengths.
  character(len=256) :: mesg    ! Message for error messages.
  real    :: dx_2 = -1.0, dy_2 = -1.0
  real    :: pi_180
  integer :: option = -1
  integer :: i, j, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq
  pi_180 = (4.0*atan(1.0))/180.0

  select case ( trim(name) )
    case ("global_1deg")    ; option = 1 ; dx_2 = 0.5*1.0
    case default ; call GOLD_error(FATAL, "reset_face_lengths_named: "//&
      "Unrecognized channel configuration name "//trim(name))
  end select

  if (option==1) then ! 1-degree settings.
    do j=jsd,jed ; do I=Isdq,Iedq  ! Change any u-face lengths within this loop.
      dy_2 = dx_2 * G%dyu(I,j)*G%Idxu(I,j) * cos(pi_180 * G%geolatu(I,j))

      if ((abs(G%geolatu(I,j)-35.5) < dy_2) .and. (G%geolonu(I,j) < -4.5) .and. &
          (G%geolonu(I,j) > -6.5)) &
        G%dy_u(I,j) = G%umask(I,j)*12000.0   ! Gibraltar

      if ((abs(G%geolatu(I,j)-12.5) < dy_2) .and. (abs(G%geolonu(I,j)-43.0) < dx_2)) &
        G%dy_u(I,j) = G%umask(I,j)*10000.0   ! Red Sea

      if ((abs(G%geolatu(i,j)-40.5) < dy_2) .and. (abs(G%geolonu(i,j)-26.0) < dx_2)) &
        G%dy_u(i,j) = G%umask(i,j)*5000.0   ! Dardanelles

      if ((abs(G%geolatu(I,j)-41.5) < dy_2) .and. (abs(G%geolonu(I,j)+220.0) < dx_2)) &
        G%dy_u(I,j) = G%umask(I,j)*35000.0   ! Tsugaru strait at 140.0e

      if ((abs(G%geolatu(I,j)-45.5) < dy_2) .and. (abs(G%geolonu(I,j)+217.5) < 0.9)) &
        G%dy_u(I,j) = G%umask(I,j)*15000.0   ! Betw Hokkaido and Sakhalin at 217&218 = 142e


      ! Greater care needs to be taken in the tripolar region.
      if ((abs(G%geolatu(I,j)-80.84) < 0.2) .and. (abs(G%geolonu(I,j)+64.9) < 0.8)) &
        G%dy_u(I,j) = G%umask(I,j)*38000.0   ! Smith Sound in Canadian Arch - tripolar region

    enddo ; enddo

    do J=Jsdq,Jedq ; do i=isd,ied  ! Change any v-face lengths within this loop.
      dy_2 = dx_2 * G%dyv(i,J)*G%Idxv(i,J) * cos(pi_180 * G%geolatv(i,J))
      if ((abs(G%geolatv(i,J)-41.0) < dy_2) .and. (abs(G%geolonv(i,J)-28.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*2500.0   ! Bosporus - should be 1000.0 m wide.

      if ((abs(G%geolatv(i,J)-13.0) < dy_2) .and. (abs(G%geolonv(i,J)-42.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*10000.0   ! Red Sea

      if ((abs(G%geolatv(i,J)+2.8) < 0.8) .and. (abs(G%geolonv(i,J)+241.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*40000.0   ! Makassar Straits at 241.5 W = 118.5 E

      if ((abs(G%geolatv(i,J)-0.56) < 0.5) .and. (abs(G%geolonv(i,J)+240.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*80000.0   ! entry to Makassar Straits at 240.5 W = 119.5 E

      if ((abs(G%geolatv(i,J)-0.19) < 0.5) .and. (abs(G%geolonv(i,J)+230.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*25000.0   ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E

      if ((abs(G%geolatv(i,J)-0.19) < 0.5) .and. (abs(G%geolonv(i,J)+229.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*25000.0   ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E

      if ((abs(G%geolatv(i,J)-0.0) < 0.25) .and. (abs(G%geolonv(i,J)+228.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*25000.0   ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E

      if ((abs(G%geolatv(i,J)+8.5) < 0.5) .and. (abs(G%geolonv(i,J)+244.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*20000.0   ! Lombok Straits at 244.5 W = 115.5 E

      if ((abs(G%geolatv(i,J)+8.5) < 0.5) .and. (abs(G%geolonv(i,J)+235.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*20000.0   ! Timor Straits at 235.5 W = 124.5 E

      if ((abs(G%geolatv(i,J)-52.5) < dy_2) .and. (abs(G%geolonv(i,J)+218.5) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*2500.0    ! Russia and Sakhalin Straits at 218.5 W = 141.5 E

      ! Greater care needs to be taken in the tripolar region.
      if ((abs(G%geolatv(i,J)-76.8) < 0.06) .and. (abs(G%geolonv(i,J)+88.7) < dx_2)) &
        G%dx_v(i,J) = G%vmask(i,J)*8400.0    ! Jones Sound in Canadian Arch - tripolar region

    enddo ; enddo
  endif

  ! These checks apply regardless of the chosen option.

  do j=jsd,jed ; do I=Isdq,Iedq
    if (G%dy_u(I,j) > G%DYu(I,j)) then
      write(mesg,'("dy_u of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
                   G%dy_u(I,j), G%DYu(I,j), G%dy_u(I,j)-G%DYu(I,j), &
                   G%geolonu(I,j), G%geolatu(I,j)
      call GOLD_error(FATAL,"reset_face_lengths_named "//mesg)
    endif
    G%dxdy_u(I,j) = G%DXu(I,j)*G%dy_u(I,j)
    G%Idxdy_u(I,j) = 0.0
    if (G%dxdy_u(I,j) > 0.0) G%Idxdy_u(I,j) = G%umask(I,j) / G%dxdy_u(I,j)
  enddo ; enddo

  do J=Jsdq,Jedq ; do i=isd,ied
    if (G%dx_v(i,J) > G%DXv(i,J)) then
      write(mesg,'("dx_v of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
                   G%dx_v(i,J), G%DXv(i,J), G%dx_v(i,J)-G%DXv(i,J), &
                   G%geolonv(i,J), G%geolatv(i,J)

      call GOLD_error(FATAL,"reset_face_lengths_named "//mesg)
    endif
    G%dxdy_v(i,J) = G%DYv(i,J)*G%dx_v(i,J)
    G%Idxdy_v(i,J) = 0.0
    if (G%dxdy_v(i,J) > 0.0) G%Idxdy_v(i,J) = G%vmask(i,J) / G%dxdy_v(i,J)
  enddo ; enddo

end subroutine reset_face_lengths_named
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------

subroutine reset_face_lengths_file(G, param_file)
  type(ocean_grid_type), intent(inout) :: G
  type(param_file_type), intent(in) :: param_file
!   This subroutine sets the open face lengths at selected points to restrict
! passages to their observed widths.

! Arguments: G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  character(len=40)  :: mod = "reset_face_lengths_file" ! This subroutine's name.
  character(len=256) :: mesg    ! Message for error messages.
  character(len=200) :: filename, chan_file, inputdir ! Strings for file/path
  integer :: i, j, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq
  ! These checks apply regardless of the chosen option.

  call GOLD_mesg("  GOLD_initialization.F90, reset_face_lengths_file: "//&
                 " setting channel configurations from file", 5)

  call get_param(param_file, mod, "CHANNEL_WIDTH_FILE", chan_file, &
                 "The file from which the list of narrowed channels is read.", &
                 default="ocean_geometry.nc")
  call get_param(param_file,  mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(chan_file)
  call log_param(param_file, mod, "INPUTDIR/CHANNEL_WIDTH_FILE", filename)

  if (is_root_pe()) then ; if (.not.file_exists(filename)) &
    call GOLD_error(FATAL," reset_face_lengths_file: Unable to open "//&
                           trim(filename))
  endif

  call read_data(filename,"dyuo",G%dy_u,domain=G%Domain%mpp_domain)
  call read_data(filename,"dxvo",G%dx_v,domain=G%Domain%mpp_domain)
  call pass_vector(G%dy_u, G%dx_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)

  do j=jsd,jed ; do I=Isdq,Iedq
    if (G%dy_u(I,j) > G%DYu(I,j)) then
      write(mesg,'("dy_u of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
                   G%dy_u(I,j), G%DYu(I,j), G%dy_u(I,j)-G%DYu(I,j), &
                   G%geolonu(I,j), G%geolatu(I,j)
      call GOLD_error(FATAL,"reset_face_lengths_file "//mesg)
    endif
    G%dxdy_u(I,j) = G%DXu(I,j)*G%dy_u(I,j)
    G%Idxdy_u(I,j) = 0.0
    if (G%dxdy_u(I,j) > 0.0) G%Idxdy_u(I,j) = G%umask(I,j) / G%dxdy_u(I,j)
  enddo ; enddo

  do J=Jsdq,Jedq ; do i=isd,ied
    if (G%dx_v(i,J) > G%DXv(i,J)) then
      write(mesg,'("dx_v of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
                   &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
                   G%dx_v(i,J), G%DXv(i,J), G%dx_v(i,J)-G%DXv(i,J), &
                   G%geolonv(i,J), G%geolatv(i,J)

      call GOLD_error(FATAL,"reset_face_lengths_file "//mesg)
    endif
    G%dxdy_v(i,J) = G%DYv(i,J)*G%dx_v(i,J)
    G%Idxdy_v(i,J) = 0.0
    if (G%dxdy_v(i,J) > 0.0) G%Idxdy_v(i,J) = G%vmask(i,J) / G%dxdy_v(i,J)
  enddo ; enddo

end subroutine reset_face_lengths_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine reset_face_lengths_list(G, param_file)
  type(ocean_grid_type), intent(inout) :: G
  type(param_file_type), intent(in) :: param_file
!   This subroutine sets the open face lengths at selected points to restrict
! passages to their observed widths.

! Arguments: G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  character(len=120), pointer, dimension(:) :: lines => NULL()
  character(len=120) :: line
  character(len=200) :: filename, chan_file, inputdir ! Strings for file/path
  character(len=40)  :: mod = "reset_face_lengths_list" ! This subroutine's name.
  real, pointer, dimension(:,:) :: &
    u_lat => NULL(), u_lon => NULL(), v_lat => NULL(), v_lon => NULL()
  real, pointer, dimension(:) :: &
    u_width => NULL(), v_width => NULL()
  real    :: lat, lon     ! The latitude and longitude of a point.
  real    :: lon_p, lon_m ! The longitude of a point shifted by 360 degrees.
  logical :: check_360    ! If true, check for longitudes that are shifted by
                          ! +/- 360 degrees from the specified range of values.
  logical :: found_u, found_v
  logical :: unit_in_use
  integer :: ios, iounit, isu, isv
  integer :: last, num_lines, nl_read, ln, npt, u_pt, v_pt
  integer :: i, j, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  call GOLD_mesg("  GOLD_initialization.F90, reset_face_lengths_list: "//&
                 " setting channel configurations from list", 5)

  call get_param(param_file, mod, "CHANNEL_LIST_FILE", chan_file, &
                 "The file from which the list of narrowed channels is read.", &
                 default="GOLD_channel_list")
  call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(chan_file)
  call log_param(param_file, mod, "INPUTDIR/CHANNEL_LIST_FILE", filename)
  call get_param(param_file, mod, "CHANNEL_LIST_360_LON_CHECK", check_360, &
                 "If true, the channel configuration list works for any \n"//&
                 "longitudes in the range of -360 to 360.", default=.true.)

  if (is_root_pe()) then
    ! Open the input file.
    if (.not.file_exists(filename)) call GOLD_error(FATAL, &
        " reset_face_lengths_list: Unable to open "//trim(filename))

    ! Find an unused unit number.
    do iounit=10,512
      INQUIRE(iounit,OPENED=unit_in_use) ; if (.not.unit_in_use) exit
    enddo
    if (iounit >= 512) call GOLD_error(FATAL, &
        "reset_face_lengths_list: No unused file unit could be found.")

    ! Open the parameter file.
    open(iounit, file=trim(filename), access='SEQUENTIAL', &
         form='FORMATTED', action='READ', position='REWIND', iostat=ios)
    if (ios /= 0) call GOLD_error(FATAL, &
            "reset_face_lengths_list: Error opening "//trim(filename))

    ! Count the number of u_width and v_width entries.
    call read_face_length_list(iounit, filename, num_lines, lines)
  endif

  ! Broadcast the number of lines and allocate the required space.
  call broadcast(num_lines, root_PE())
  u_pt = 0 ; v_pt = 0
  if (num_lines > 0) then
    allocate (lines(num_lines))
    if (num_lines > 0) then
      allocate(u_lat(2,num_lines)) ; u_lat(:,:) = -1e34
      allocate(u_lon(2,num_lines)) ; u_lon(:,:) = -1e34
      allocate(u_width(num_lines)) ; u_width(:) = -1e34

      allocate(v_lat(2,num_lines)) ; v_lat(:,:) = -1e34
      allocate(v_lon(2,num_lines)) ; v_lon(:,:) = -1e34
      allocate(v_width(num_lines)) ; v_width(:) = -1e34
    endif

    ! Actually read the lines.
    if (is_root_pe()) then
      call read_face_length_list(iounit, filename, nl_read, lines)
      if (nl_read /= num_lines) &
        call GOLD_error(FATAL, 'reset_face_lengths_list : Found different '// &
                  'number of valid lines on second reading of '//trim(filename))
      close(iounit) ; iounit = -1
    endif

    ! Broadcast the lines.
    call broadcast(lines, 120, root_PE())

    ! Populate the u_width, etc., data.
    do ln=1,num_lines
      line = lines(ln)
      ! Detect keywords
      found_u = .false.; found_v = .false.
      isu = index(uppercase(line), "U_WIDTH" ); if (isu > 0) found_u = .true.
      isv = index(uppercase(line), "V_WIDTH" ); if (isv > 0) found_v = .true.
      
      ! Store and check the relevant values.
      if (found_u) then
        u_pt = u_pt + 1
        read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt) 
        if (is_root_PE()) then
          if (check_360) then
            if ((abs(u_lon(1,u_pt)) > 360.0) .or. (abs(u_lon(2,u_pt)) > 360.0)) &
              call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-bounds "//&
                 "u-longitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
            if ((abs(u_lat(1,u_pt)) > 180.0) .or. (abs(u_lat(2,u_pt)) > 180.0)) &
              call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-bounds "//&
                 "u-latitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
          endif
          if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
            call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-order "//&
               "u-face latitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
            call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-order "//&
               "u-face longitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (u_width(u_pt) < 0.0) &
            call GOLD_error(WARNING, "reset_face_lengths_list : Negative "//&
               "u-width found when reading line "//trim(line)//" from file "//&
               trim(filename))
        endif
      elseif (found_v) then
        v_pt = v_pt + 1
        read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
        if (is_root_PE()) then
          if (check_360) then
            if ((abs(v_lon(1,v_pt)) > 360.0) .or. (abs(v_lon(2,v_pt)) > 360.0)) &
              call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-bounds "//&
                 "v-longitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
            if ((abs(v_lat(1,v_pt)) > 180.0) .or. (abs(v_lat(2,v_pt)) > 180.0)) &
              call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-bounds "//&
                 "v-latitude found when reading line "//trim(line)//" from file "//&
                 trim(filename))
          endif
          if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
            call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-order "//&
               "v-face latitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
            call GOLD_error(WARNING, "reset_face_lengths_list : Out-of-order "//&
               "v-face longitudes found when reading line "//trim(line)//" from file "//&
               trim(filename))
          if (v_width(v_pt) < 0.0) &
            call GOLD_error(WARNING, "reset_face_lengths_list : Negative "//&
               "v-width found when reading line "//trim(line)//" from file "//&
               trim(filename))
        endif
      endif
    enddo
    
    deallocate(lines)
  endif

  do j=jsd,jed ; do I=Isdq,Iedq
    lat = G%geolatu(I,j) ; lon = G%geolonu(I,j)
    if (check_360) then ; lon_p = lon+360.0 ; lon_m = lon-360.0
    else ; lon_p = lon ; lon_m = lon ; endif

    do npt=1,u_pt
      if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
          (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
           ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
           ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) &
        G%dy_u(I,j) = G%umask(I,j) * min(G%DYu(I,j), max(u_width(npt), 0.0))
    enddo

    G%dxdy_u(I,j) = G%DXu(I,j)*G%dy_u(I,j)
    G%Idxdy_u(I,j) = 0.0
    if (G%dxdy_u(I,j) > 0.0) G%Idxdy_u(I,j) = G%umask(I,j) / G%dxdy_u(I,j)
  enddo ; enddo

  do J=Jsdq,Jedq ; do i=isd,ied
    lat = G%geolatv(i,J) ; lon = G%geolonv(i,J)
    if (check_360) then ; lon_p = lon+360.0 ; lon_m = lon-360.0
    else ; lon_p = lon ; lon_m = lon ; endif

    do npt=1,v_pt
      if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
          (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
           ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
           ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) &
        G%dx_v(i,J) = G%vmask(i,J) * min(G%DXv(i,J), max(v_width(npt), 0.0))
    enddo

    G%dxdy_v(i,J) = G%DYv(i,J)*G%dx_v(i,J)
    G%Idxdy_v(i,J) = 0.0
    if (G%dxdy_v(i,J) > 0.0) G%Idxdy_v(i,J) = G%vmask(i,J) / G%dxdy_v(i,J)
  enddo ; enddo

  if (num_lines > 0) then
    deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width)
    deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width)
  endif

end subroutine reset_face_lengths_list

subroutine read_face_length_list(iounit, filename, num_lines, lines)
  integer,               intent(in)  :: iounit
  character(len=*),      intent(in)  :: filename
  integer,               intent(out) :: num_lines
  character(len=120), pointer, dimension(:) :: lines

  !   This subroutine reads and counts the non-blank lines in the face length
  ! list file, after removing comments.
  character(len=120) :: line, line_up
  logical :: found_u, found_v
  integer :: isu, isv, icom, verbose
  integer :: last
  
  num_lines = 0

  if (iounit <= 0) return
  rewind(iounit)
  do while(.true.)
    read(iounit, '(a)', end=8, err=9) line
    last = len_trim(line)
    ! Eliminate either F90 or C comments from the line.
    icom = index(line(:last), "!") ; if (icom > 0) last = icom-1
    icom = index(line(:last), "/*") ; if (icom > 0) last = icom-1
    if (last < 1) cycle

    ! Detect keywords
    line_up = uppercase(line)
    found_u = .false.; found_v = .false.
    isu = index(line_up(:last), "U_WIDTH" ); if (isu > 0) found_u = .true.
    isv = index(line_up(:last), "V_WIDTH" ); if (isv > 0) found_v = .true.

    if (found_u .and. found_v) call GOLD_error(FATAL, &
      "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
      "reading the line "//trim(line(:last))//" in file "//trim(filename))
    if (found_u .or. found_v) then
      num_lines = num_lines + 1
      if (associated(lines)) then
        lines(num_lines) = line(1:last)
      endif
    endif
  enddo ! while (.true.)

8 continue
  return

9 call GOLD_error(FATAL, "read_face_length_list : "//&
                  "Error while reading file "//trim(filename))
 
end subroutine read_face_length_list

! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_velocity_depth_max(G)
  type(ocean_grid_type), intent(inout) :: G
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! maximum of the adjacent depths.
  integer :: i, j

  do I=G%isd,G%ied-1 ; do j=G%jsd,G%jed
    G%Dblock_u(I,j) = G%umask(I,j)*max(G%D(i,j),G%D(i+1,j))
    G%Dopen_u(I,j) = G%Dblock_u(I,j)
  enddo ; enddo
  do i=G%isd,G%ied ; do J=G%jsd,G%jed-1
    G%Dblock_v(I,J) = G%vmask(i,J)*max(G%D(i,j),G%D(i,j+1))
    G%Dopen_v(I,J) = G%Dblock_v(I,J)
  enddo ; enddo
end subroutine set_velocity_depth_max
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine set_velocity_depth_min(G)
  type(ocean_grid_type), intent(inout) :: G
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! minimum of the adjacent depths.
  integer :: i, j

  do I=G%isd,G%ied-1 ; do j=G%jsd,G%jed
    G%Dblock_u(I,j) = G%umask(I,j)*min(G%D(i,j),G%D(i+1,j))
    G%Dopen_u(I,j) = G%Dblock_u(I,j)
  enddo ; enddo
  do i=G%isd,G%ied ; do J=G%jsd,G%jed-1
    G%Dblock_v(I,J) = G%vmask(i,J)*min(G%D(i,j),G%D(i,j+1))
    G%Dopen_v(I,J) = G%Dblock_v(I,J)
  enddo ; enddo
end subroutine set_velocity_depth_min
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine write_ocean_geometry_file(G, param_file, directory)
  type(ocean_grid_type), intent(inout) :: G
  type(param_file_type), intent(in) :: param_file
  character(len=*),      intent(in) :: directory
!   This subroutine writes out a file containing all of the ocean geometry
! and grid data uses by the GOLD ocean model.
! Arguments: G - The ocean's grid structure.  Effectively intent in.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      directory - The directory into which to place the file.
  character(len=120) :: filepath
  character(len=40)  :: mod = "write_ocean_geometry_file"
  integer, parameter :: nFlds=23
  type(vardesc) :: vars(nFlds)
  type(fieldtype) :: fields(nFlds)
  integer :: unit
  integer :: file_threading
  integer :: nFlds_used
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
  logical :: multiple_files
  real :: out_h(SZI_(G),SZJ_(G))
  real :: out_u(SZIQ_(G),SZJ_(G))
  real :: out_v(SZI_(G),SZJQ_(G))
  real :: out_q(SZIQ_(G),SZJQ_(G))
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

!   vardesc is a structure defined in GOLD_io.F90.  The elements of
! this structure, in order, are:
! (1) the variable name for the NetCDF file
! (2) the variable's long name
! (3) a character indicating the  horizontal grid, which may be '1' (column),
!     'h', 'q', 'u', or 'v', for the corresponding C-grid variable
! (4) a character indicating the vertical grid, which may be 'L' (layer),
!     'i' (interface), or '1' (no vertical location)
! (5) a character indicating the time levels of the field, which may be
!    's' (snap-shot), 'a' (average between snapshots), 'm' (monthly average),
!     or '1' (no time variation) between snapshots
! (6) the variable's units
! (7) a character indicating the size in memory to write, which may be
!     'd' (8-byte) or 'f' (4-byte).
  vars(1) = vardesc("geolatb","latitude at q points",'q','1','1',"degree",'d')
  vars(2) = vardesc("geolonb","longitude at q points",'q','1','1',"degree",'d')
  vars(3) = vardesc("geolat", "latitude at h points", 'h','1','1',"degree",'d')
  vars(4) = vardesc("geolon","longitude at h points",'h','1','1',"degree",'d')
  vars(5) = vardesc("D","Basin Depth",'h','1','1',"meter",'d')
  vars(6) = vardesc("f","Coriolis Parameter",'q','1','1',"second-1", 'd')
  vars(7) = vardesc("dxv","Zonal grid spacing at v points",'v','1','1',"m",'d')
  vars(8) = vardesc("dyu","Meridional grid spacing at u points",'u','1','1',"m",'d')
  vars(9) = vardesc("dxu","Zonal grid spacing at u points",'u','1','1',"m",'d')
  vars(10)= vardesc("dyv","Meridional grid spacing at v points",'v','1','1',"m",'d')
  vars(11)= vardesc("dxh","Zonal grid spacing at h points",'h','1','1',"m",'d')
  vars(12)= vardesc("dyh","Meridional grid spacing at h points",'h','1','1',"m",'d')
  vars(13)= vardesc("dxq","Zonal grid spacing at q points",'q','1','1',"m",'d')
  vars(14)= vardesc("dyq","Meridional grid spacing at q points",'q','1','1',"m",'d')
  vars(15)= vardesc("Ah","Area of h cells",'h','1','1',"m2",'d')
  vars(16)= vardesc("Aq","Area of q cells",'q','1','1',"m2",'d')

  vars(17)= vardesc("dxvo","Open zonal grid spacing at v points",'v','1','1',"m",'d')
  vars(18)= vardesc("dyuo","Open meridional grid spacing at u points",'u','1','1',"m",'d')
  vars(19)= vardesc("wet", "land or ocean?", 'h','1','1',"none",'d')

  vars(20) = vardesc("Dblock_u","Blocked depth at u points",'u','1','1',"meter",'d')
  vars(21) = vardesc("Dopen_u","Open depth at u points",'u','1','1',"meter",'d')
  vars(22) = vardesc("Dblock_v","Blocked depth at v points",'v','1','1',"meter",'d')
  vars(23) = vardesc("Dopen_v","Open depth at v points",'v','1','1',"meter",'d')

  nFlds_used = 19 ; if (G%bathymetry_at_vel) nFlds_used = 23

  filepath = trim(directory) // "ocean_geometry"

  out_h(:,:) = 0.0
  out_u(:,:) = 0.0
  out_v(:,:) = 0.0
  out_q(:,:) = 0.0

  call get_param(param_file, mod, "PARALLEL_RESTARTFILES", multiple_files, &
                 "If true, each processor writes its own restart file, \n"//&
                 "otherwise a single restart file is generated", &
                 default=.false.)
  file_threading = SINGLE_FILE
  if (multiple_files) file_threading = MULTIPLE

  call create_file(unit, trim(filepath), vars, nFlds_used, G, fields, file_threading)

  do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geolatq(I,J); enddo; enddo
  call write_field(unit, fields(1), G%Domain%mpp_domain, out_q)
  do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geolonq(I,J); enddo; enddo
  call write_field(unit, fields(2), G%Domain%mpp_domain, out_q)
  call write_field(unit, fields(3), G%Domain%mpp_domain, G%geolath)
  call write_field(unit, fields(4), G%Domain%mpp_domain, G%geolonh)

  call write_field(unit, fields(5), G%Domain%mpp_domain, G%D)
  call write_field(unit, fields(6), G%Domain%mpp_domain, G%f)

  do J=Jsq,Jeq; do i=is,ie; out_v(i,J) = G%DXv(i,J); enddo; enddo
  call write_field(unit, fields(7), G%Domain%mpp_domain, out_v)
  do j=js,je; do I=Isq,Ieq; out_u(I,j) = G%DYu(I,j); enddo; enddo
  call write_field(unit, fields(8), G%Domain%mpp_domain, out_u)

  do J=Jsq,Jeq; do i=is,ie; out_u(i,J) = G%DXu(i,J); enddo; enddo
  call write_field(unit, fields(9), G%Domain%mpp_domain, out_u)
  do j=js,je; do I=Isq,Ieq; out_v(I,j) = G%DYv(I,j); enddo; enddo
  call write_field(unit, fields(10), G%Domain%mpp_domain, out_v)

  do J=Jsq,Jeq; do i=is,ie; out_h(i,J) = G%DXh(i,J); enddo; enddo
  call write_field(unit, fields(11), G%Domain%mpp_domain, out_h)
  do j=js,je; do I=Isq,Ieq; out_h(I,j) = G%DYh(I,j); enddo; enddo
  call write_field(unit, fields(12), G%Domain%mpp_domain, out_h)

  do J=Jsq,Jeq; do i=is,ie; out_q(i,J) = G%DXq(i,J); enddo; enddo
  call write_field(unit, fields(13), G%Domain%mpp_domain, out_q)
  do j=js,je; do I=Isq,Ieq; out_q(I,j) = G%DYq(I,j); enddo; enddo
  call write_field(unit, fields(14), G%Domain%mpp_domain, out_q)

  do j=js,je; do i=is,ie; out_h(i,j) = G%DXDYh(i,j); enddo; enddo
  call write_field(unit, fields(15), G%Domain%mpp_domain, out_h)
  do j=js,je; do i=is,ie; out_q(i,j) = G%DXDYq(i,j); enddo; enddo
  call write_field(unit, fields(16), G%Domain%mpp_domain, out_q)

!  do J=Jsq,Jeq; do i=is,ie; out_v(i,J) = G%dx_v(i,J); enddo; enddo
  call write_field(unit, fields(17), G%Domain%mpp_domain, G%dx_v)
!  do j=js,je; do I=Isq,Ieq; out_u(I,j) = G%dy_u(I,j); enddo; enddo
  call write_field(unit, fields(18), G%Domain%mpp_domain, G%dy_u)
  call write_field(unit, fields(19), G%Domain%mpp_domain, G%hmask)

  if (G%bathymetry_at_vel) then
    call write_field(unit, fields(20), G%Domain%mpp_domain, G%Dblock_u)
    call write_field(unit, fields(21), G%Domain%mpp_domain, G%Dopen_u)
    call write_field(unit, fields(22), G%Domain%mpp_domain, G%Dblock_v)
    call write_field(unit, fields(23), G%Domain%mpp_domain, G%Dopen_v)
  endif

  call close_file(unit)

end subroutine write_ocean_geometry_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine write_vertgrid_file(G, param_file, directory)
  type(ocean_grid_type), intent(inout) :: G
  type(param_file_type), intent(in) :: param_file
  character(len=*),      intent(in) :: directory
!   This subroutine writes out a file containing any available data related
! to the vertical grid used by the GOLD ocean model.
! Arguments: G - The ocean's grid structure.  Effectively intent in.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      directory - The directory into which to place the file.
  character(len=120) :: filepath
  type(vardesc) :: vars(2)
  type(fieldtype) :: fields(2)
  integer :: unit

  filepath = trim(directory) // trim("Vertical_coordinate")

  vars(1) = vardesc("R","Target Potential Density",'1','L','1',"kilogram meter-3", 'd')
  vars(2) = vardesc("g","Reduced gravity",'1','L','1',"meter second-2", 'd')

  call create_file(unit, trim(filepath), vars, 2, G, fields, SINGLE_FILE)

  call write_field(unit, fields(1), G%Rlay)
  call write_field(unit, fields(2), G%g_prime)

  call close_file(unit)

end subroutine write_vertgrid_file
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
subroutine Get_GOLD_Input(param_file, dirs, check_params)
  use GOLD_io, only : open_namelist_file, check_nml_error
  type(param_file_type), optional, intent(out) :: param_file
  type(directories),     optional, intent(out) :: dirs
  logical,               optional, intent(in)  :: check_params

!    See if the run is to be started from saved conditions, and get  !
!  the names of the I/O directories and initialization file.  This   !
!  subroutine also calls the subroutine that allows run-time changes !
!  in parameters.                                                    !
!    Read from a namelist instead of getting strings from the command line as
!  in the C-version of GOLD.
  integer, parameter :: npf = 5 ! Maximum number of parameter files
  character(len=120) :: &
    parameter_filename(npf) = ' ', & ! List of files containing parameters.
    output_directory = ' ', &   ! Directory to use to write the model output.
    restart_input_dir = ' ', &  ! Directory for reading restart and input files.
    restart_output_dir = ' ', & ! Directory into which to write restart files.
    input_filename  = ' '       ! A string that indicates the input files or how
                                ! the run segment should be started.
  integer :: unit, io, ierr

  namelist /GOLD_input_nml/ output_directory, input_filename, parameter_filename, &
                           restart_input_dir, restart_output_dir

  if (file_exists('input.nml')) then
    unit = open_namelist_file(file='input.nml')
  else
    call GOLD_error(FATAL,'Required namelist file input.nml does not exist.')
  endif

  ierr=1 ; do while (ierr /= 0)
    read(unit, nml=GOLD_input_nml, iostat=io, end=10)
    ierr = check_nml_error(io, 'GOLD_input_nml')
  enddo
10 call close_file(unit)
  if (present(dirs)) then
    dirs%output_directory = slasher(output_directory)
    dirs%restart_output_dir = slasher(restart_output_dir)
    dirs%restart_input_dir = slasher(restart_input_dir)
    dirs%input_filename = input_filename
  endif

  if (present(param_file)) then ; do io = 1, npf
    if (len_trim(trim(parameter_filename(io))) > 0) &
      call open_param_file(trim(parameter_filename(io)), param_file, check_params)
  enddo ; endif

end subroutine Get_GOLD_Input
! -----------------------------------------------------------------------------

subroutine GOLD_temp_salt_initialize_from_Z(h, tv, G, PF, dirs)
  real, dimension(NXMEM_,NYMEM_,NZ_), intent(out)   :: h    
  type(thermo_var_ptrs),              intent(inout) :: tv
  type(ocean_grid_type),              intent(in)    :: G
  type(param_file_type),              intent(in)    :: PF
  type(directories),                  intent(in)    :: dirs
!   This subroutine determines the isopycnal interfaces and layer potential
! temperatures and salinities directly from a z-space file on a latitude-
! longitude grid.
!
!   This subroutine was written by M. Harrison, with input from R. Hallberg.
!
! Arguments: 
!  (out)     h  - Layer thickness, in m.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (inout)   G       - The ocean's grid structure.
!  (in)      PF      - A structure indicating the open file to parse for
!                      model parameter values.
!  (in)      dirs    - A structure containing several relevant directory paths.


  character(len=200) :: filename   ! The name of an input file containing temperature
                                   ! and salinity in z-space.
  character(len=200) :: inputdir ! The directory where NetCDF input files are.
  character(len=200) :: mesg

  type(EOS_type), pointer :: eos => NULL()
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod = "GOLD_initialize_layers_from_Z" ! This module's name.

  integer :: is, ie, js, je, nz ! compute domain indices
  integer :: isc,iec,jsc,jec    ! global compute domain indices
  integer :: isg, ieg, jsg, jeg ! global extent
  integer :: isd, ied, jsd, jed ! data domain indices
  integer :: rcode, no_fill
  integer :: ndims                 

  integer :: nkml, nkbl         ! number of mixed and buffer layers
  integer :: xhalo,yhalo        ! halo widths
  real    :: PI_180             ! for conversion from degrees to radians

  type(horiz_interp_type) :: Interp
  logical :: new_sim
  logical :: correct_thickness

  character(len=40) :: potemp_var, salin_var
  character(len=8)  :: laynum
  integer, parameter :: niter=10   ! number of iterations for t/s adjustment to layer density

  real, parameter    :: missing_value=-1.e34
  real, parameter    :: temp_land_fill = 0.0, salt_land_fill=35.0

  !data arrays
  real, dimension(:,:), allocatable :: x_in, y_in, nlevs
  real, dimension(:), allocatable :: lon_in, lon_in_p, lat_in, lat_in_p, last_row
  real, dimension(:), allocatable :: z_in, z_edges_in, Rb
  real, dimension(:,:), allocatable :: temp_in, salt_in, mask_in
  real :: npole, pole, max_lat, min_depth, dilate
  integer :: ncid, id, jd, kd, varid_t, varid_s, jdp, inconsistent
  integer, dimension(4) :: start, count, dims    
  real :: missing_value_temp, missing_value_salt    
  real, dimension(:,:), allocatable :: tmp_in ! A 2-d array for holding input data.  
  !global arrays 
  real, dimension(:,:), allocatable :: temp_prev, salt_prev, mask_prev    
  real, dimension(:,:), allocatable :: temp_out, salt_out, rho_out, mask_out
  real, dimension(:,:), allocatable :: hmask, Depth
  real, dimension(:,:), allocatable :: lon_out, lat_out, x_out, y_out
  integer, dimension(:,:), allocatable :: good, fill
  real, dimension(:), allocatable   :: press
  real, dimension(:,:), allocatable :: hlay

  !local arrays
  real, dimension(:,:,:), allocatable :: temp_z, salt_z, mask_z, rho_z
  real, dimension(:,:,:), allocatable :: zi

  integer ::  i, j, k, ks, np, ni, nj

  logical :: reentrant_x, tripolar_n, add_np
  logical :: debug_point = .false.  ! manually set this to true and adjust the locations below
                                    ! if you want to track down a problem with a particular
                                    ! location.
  logical, parameter :: adjust_temperature = .true.

  real, parameter :: debug_lon=-180.5, debug_lat=0.125
  real, dimension(:,:,:), allocatable :: zip    

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg

  PI_180=atan(1.0)/45.

  call GOLD_mesg(" GOLD_temp_salt_initiale_from_Z.F90, GOLD_temp_salt_initialize_from_Z: subroutine entered", 3)
  call log_version(PF, mod, version, tagname)

  new_sim = .false.
  if ((dirs%input_filename(1:1) == 'n') .and. &
       (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true.

  inputdir = "." ;  call get_param(PF, mod, "INPUTDIR", inputdir)
  inputdir = slasher(inputdir)    

  eos => tv%eqn_of_state

  call mpp_get_compute_domain(G%domain%mpp_domain,isc,iec,jsc,jec)

  reentrant_x = .false. ;  call get_param(PF, mod, "REENTRANT_X", reentrant_x,default=.true.)
  tripolar_n = .false. ;  call get_param(PF, mod, "TRIPOLAR_N", tripolar_n, default=.false.)
  call get_param(PF, mod, "MINIMUM_DEPTH", min_depth, default=0.0)

  call get_param(PF, mod, "NKML",nkml)
  call get_param(PF, mod, "NKBL",nkbl)    

  call get_param(PF, mod, "TEMP_SALT_Z_INIT_FILE",filename, &
                 "The name of the z-space input file used to initialize \n"//&
                 "the layer thicknesses, temperatures and salinities.", &
                 default="temp_salt_z.nc")
  filename = trim(inputdir)//trim(filename)
  call get_param(PF, mod, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
                 "The name of the potential temperature variable in \n"//&
                 "TEMP_SALT_Z_INIT_FILE.", default="ptemp")
  call get_param(PF, mod, "Z_INIT_FILE_SALT_VAR", salin_var, &
                 "The name of the salinity variable in \n"//&
                 "TEMP_SALT_Z_INIT_FILE.", default="salt")

!   Read input grid coordinates for temperature and salinity field
!   in z-coordinate dataset. The file is REQUIRED to contain the
!   following:
!
!   dimension variables:
!            lon (degrees_E), lat (degrees_N), depth(meters)
!   variables:
!            ptemp(lon,lat,depth) : degC, potential temperature
!            salt (lon,lat,depth) : PSU, salinity
!
!   The first record will be read if there are multiple time levels.
!   The observation grid MUST tile the model grid. If the model grid extends
!   to the North Pole, the input data are extrapolated using the average
!   value at the northernmost latitude.      

  rcode = NF90_OPEN(filename, NF90_NOWRITE, ncid)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error opening file "//trim(filename)//&
                           " in GOLD_initialize_layers_from_Z")
  rcode = NF90_INQ_VARID(ncid, potemp_var, varid_t)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error finding variable "//trim(potemp_var)//&
                                 " in file "//trim(filename)//" in GOLD_initialize_layers_from_Z")        
  rcode = NF90_INQ_VARID(ncid, salin_var, varid_s)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error finding variable "//trim(salin_var)//&
                                 " in file "//trim(filename)//" in GOLD_initialize_layers_from_Z")        
  rcode = NF90_INQUIRE_VARIABLE(ncid, varid_t, ndims=ndims, dimids=dims)
  if (rcode .ne. 0) call GOLD_error(FATAL,'error inquiring dimensions GOLD_initialize_layers_from_Z')            
  if (ndims < 3) call GOLD_error(FATAL,"Variable "//trim(potemp_var)//" in file "// &
              trim(filename)//" has too few dimensions.")            
  rcode = NF90_INQUIRE_DIMENSION(ncid, dims(1), len=id)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error reading dimension 1 data for "// &
                trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")               
  rcode = NF90_INQUIRE_DIMENSION(ncid, dims(2), len=jd)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error reading dimension 2 data for "// &
                trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")               
  rcode = NF90_INQUIRE_DIMENSION(ncid, dims(3), len=kd)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error reading dimension 3 data for "// &
                trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")               

  missing_value_temp=0.0 ; missing_value_salt=0.0
  rcode = NF90_GET_ATT(ncid, varid_t, "_FillValue", missing_value_temp)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error finding missing value for "//&
       trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")    
  rcode = NF90_GET_ATT(ncid, varid_s, "_FillValue", missing_value_salt)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error finding missing value for "//&
       trim(salin_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")    

  allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1))
  allocate(temp_z(is:ie,js:je,kd), salt_z(is:ie,js:je,kd), rho_z(is:ie,js:je,kd), mask_z(is:ie,js:je,kd))

  start = 1; count = 1; count(1) = id
  rcode = NF90_GET_VAR(ncid, dims(1), lon_in, start, count)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error reading dimension 1 values for "// &
                trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")               
  start = 1; count = 1; count(1) = jd
  rcode = NF90_GET_VAR(ncid, dims(2), lat_in, start, count)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error reading dimension 2 values for "// &
                trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")               
  start = 1; count = 1; count(1) = kd
  rcode = NF90_GET_VAR(ncid, dims(3), z_in, start, count)
  if (rcode .ne. 0) call GOLD_error(FATAL,"error reading dimension 3 values for "// &
                trim(potemp_var)//" in file "// trim(filename)//" in GOLD_initialize_layers_from_Z")               

! extrapolate the input data to the north pole using the northerm-most latitude

  max_lat = maxval(lat_in)
  add_np=.false.
  if (max_lat < 90.0) then
    add_np=.true.
    jdp=jd+1
    allocate(lat_in_p(jdp))
    lat_in_p(1:jd)=lat_in(:)
    lat_in_p(jd+1)=90.0
    deallocate(lat_in)
    allocate(lat_in(1:jdp))
    lat_in(:)=lat_in_p(:)
  else
    jdp=jd
  endif

! construct level cell boundaries as the mid-point between adjacent centers

  do k=2,kd
   z_edges_in(k)=0.5*(z_in(k-1)+z_in(k))
  enddo
  z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1)

  call horiz_interp_init()

  lon_in = lon_in*PI_180
  lat_in = lat_in*PI_180
  allocate(x_in(id,jdp),y_in(id,jdp))        
  call meshgrid(lon_in,lat_in, x_in, y_in)

  allocate(lon_out(isd:ied,jsd:jed),lat_out(isd:ied,jsd:jed))
  lon_out(is:ie,js:je) = G%geolonh(is:ie,js:je)*PI_180
  lat_out(is:ie,js:je) = G%geolath(is:ie,js:je)*PI_180

! get the global model grid
  ni=ieg-isg+1 ; nj = jeg-jsg+1
  allocate(x_out(ni,nj),y_out(ni,nj))
  call mpp_global_field(G%domain%mpp_domain,lon_out,x_out)
  call mpp_global_field(G%domain%mpp_domain,lat_out,y_out)    


  allocate(tmp_in(id,jd)) ; tmp_in(:,:)=0.0
  allocate(temp_in(id,jdp)) ; temp_in(:,:)=0.0
  allocate(salt_in(id,jdp)) ; salt_in(:,:)=0.0
  allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0
  allocate(last_row(id))    ; last_row(:)=0.0
  allocate(temp_out(ni,nj),salt_out(ni,nj),mask_out(ni,nj))
  allocate(temp_prev(ni,nj),salt_prev(ni,nj),rho_out(ni,nj))
  allocate(fill(ni,nj),hmask(ni,nj),good(ni,nj), Depth(ni,nj))
  allocate(press(ni)) ; press(:)=tv%p_ref

  temp_prev=0.0; salt_prev=0.0

! get the global wet mask and depth arrays
  call mpp_global_field(G%domain%mpp_domain,G%hmask,hmask)
  call mpp_global_field(G%domain%mpp_domain,G%D,Depth)    


! loop through each data level and interpolate to model grid.
! after interpolating, fill in points which will be needed
! to define the layers

  do k=1,kd
    write(laynum,'(I8)') k ; laynum = adjustl(laynum)
  
    start = 1; start(3) = k; count = 1; count(1) = id; count(2) = jd
    rcode = NF90_GET_VAR(ncid,varid_t, tmp_in, start, count)
    if (rcode .ne. 0) call GOLD_error(FATAL,"GOLD_initialize_layers_from_Z: "//&
                "error reading level "//trim(laynum)//" of variable "//&
                trim(potemp_var)//" in file "// trim(filename))

    if (add_np) then
      last_row(:)=tmp_in(:,jd); pole=0.0;npole=0.0
      do i=1,id
        if (abs(tmp_in(i,jd)-missing_value_temp) .gt. abs(G%Angstrom_Z*missing_value_temp)) then
          pole = pole+last_row(i)
          npole = npole+1.0
        endif
      enddo
      if (npole > 0) then
        pole=pole/npole
      else
        pole=missing_value
      endif
      temp_in(:,1:jd) = tmp_in(:,:)
      temp_in(:,jdp) = pole
    else
      temp_in(:,:) = tmp_in(:,:)
    endif

    rcode = NF90_GET_VAR(ncid, varid_s, tmp_in, start, count)
    if (rcode .ne. 0) call GOLD_error(FATAL,"GOLD_initialize_layers_from_Z: "//&
                "error reading level "//trim(laynum)//" of variable "//&
                trim(salin_var)//" in file "// trim(filename))

    if (add_np) then
      last_row(:)=tmp_in(:,jd); pole=0.0;npole=0.0
      do i=1,id
        if (abs(tmp_in(i,jd)-missing_value_salt) .gt. abs(G%Angstrom_Z*missing_value_salt)) then              
          pole = pole+last_row(i)
          npole = npole+1.0
        endif
      enddo
      if (npole > 0) then
        pole = pole/npole
      else
        pole = missing_value
      endif           
      salt_in(:,1:jd) = tmp_in(:,:)
      salt_in(:,jdp) = pole
    else
      salt_in(:,:) = tmp_in(:,:)
    endif

    mask_in=0.0

    do j=1,jdp
      do i=1,id
        if (abs(temp_in(i,j)-missing_value_temp) .gt. abs(G%Angstrom_Z*missing_value_temp)) then                           
            mask_in(i,j)=1.0
        endif
      enddo
    enddo

! call fms routine horiz_interp to interpolate input level data to model horizontal grid

    if (k .eq. 1) then
      call horiz_interp_new(Interp,x_in,y_in,x_out,y_out, &
               interp_method='bilinear',src_modulo=reentrant_x, &
               mask_in=mask_in,mask_out=mask_out)
    endif

    call horiz_interp(Interp,temp_in,temp_out,mask_in=mask_in, &
                      mask_out=mask_out,missing_value=missing_value)
    call horiz_interp(Interp,salt_in,salt_out,mask_in=mask_in, &
                      mask_out=mask_out,missing_value=missing_value)       

    fill = 0; good = 0

    do j=1,nj
      do i=1,ni
        if (mask_out(i,j).lt.1.0) then
          temp_out(i,j)=missing_value
          salt_out(i,j)=missing_value
        else
          good(i,j)=1
        endif
        if (hmask(i,j) .eq. 1.0 .and. z_edges_in(k) <= Depth(i,j) .and. mask_out(i,j) .lt. 1.0) fill(i,j)=1
      enddo
    enddo


! temp_out,salt_out contain input z-space data on the model grid with missing values
! now fill in missing values using "ICE-nine" algorithm. 

    temp_out = fill_miss_2d(temp_out,good,fill,temp_prev,reentrant_x,tripolar_n)
    salt_out = fill_miss_2d(salt_out,good,fill,salt_prev,reentrant_x,tripolar_n)

    good=good+fill

    temp_out=temp_out*hmask
    salt_out=salt_out*hmask

    temp_prev=temp_out
    salt_prev=salt_out

! next use the equation of state to create a potential density field using filled z-data

    do j=1,nj
      call calculate_density(temp_out(:,j),salt_out(:,j), press, rho_out(:,j), 1, ni, eos)
    enddo

! copy to local arrays
    xhalo=is-isd;yhalo=js-jsd
    temp_z(is:ie,js:je,k) = temp_out(isc-xhalo:iec-xhalo,jsc-yhalo:jec-yhalo)
    salt_z(is:ie,js:je,k) = salt_out(isc-xhalo:iec-xhalo,jsc-yhalo:jec-yhalo)
    mask_z(is:ie,js:je,k) = good(isc-xhalo:iec-xhalo,jsc-yhalo:jec-yhalo)
    rho_z(is:ie,js:je,k) = rho_out(isc-xhalo:iec-xhalo,jsc-yhalo:jec-yhalo)       

  enddo

  call horiz_interp_del(Interp)    

! Done with horizontal interpolation.    
! next find interface positions using local arrays
! nlevs contains the number of valid data points in each column

  allocate(zi(is:ie,js:je,nz+1))
  allocate(nlevs(is:ie,js:je))
  nlevs = sum(mask_z(is:ie,js:je,:),dim=3)

! Rb contains the layer interface densities

  allocate(Rb(nz+1))
  do k=2,nz
     Rb(k)=0.5*(G%Rlay(k-1)+G%Rlay(k))
  enddo
  Rb(1)=0.0
  Rb(nz+1)=2.0*G%Rlay(nz) - G%Rlay(nz-1)

  zi = find_interfaces(rho_z(is:ie,js:je,:),z_in,Rb,G%D(is:ie,js:je),nlevs,nkml,nkbl,min_depth)

  call get_param(PF,mod,"ADJUST_THICKNESS",correct_thickness,default=.false.)

  if (correct_thickness) then
    ! All mass below the bottom removed if the topography is shallower than
    ! the input file would indicate.  G%D is positive downward,
    ! eta is negative downward.
    do j=js,je ; do i=is,ie
      if (-zi(i,j,nz+1) > G%D(i,j) + 0.1) zi(i,j,nz+1) = -G%D(i,j)
    enddo ; enddo
  endif

  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
    if (zi(i,j,K) < (zi(i,j,K+1) + G%Angstrom_z)) then
      zi(i,j,K) = zi(i,j,K+1) + G%Angstrom_z
      h(i,j,k) = G%Angstrom_z
    else
      h(i,j,k) = zi(i,j,K) - zi(i,j,K+1)
    endif
  enddo ; enddo ; enddo


!  Check for consistency between the interface heights and topography.!
  if (correct_thickness) then
    do j=js,je ; do i=is,ie
      !   The whole column is dilated to accomodate deeper topography than
      ! the input file would indicate.
      if (-zi(i,j,nz+1) < G%D(i,j) - 0.1) then
        dilate = (zi(i,j,1)+G%D(i,j)) / (zi(i,j,1)-zi(i,j,nz+1))
        do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
      endif
    enddo ; enddo
  else
    inconsistent=0
    do j=js,je ; do i=is,ie
      if (abs(zi(i,j,nz+1) + G%D(i,j)) > 1.0) inconsistent = inconsistent + 1
    enddo ; enddo
    call sum_across_PEs(inconsistent)

    if ((inconsistent > 0) .and. (is_root_pe())) then
      write(mesg,'("Thickness initial conditions are inconsistent ",'// &
               '"with topography in ",I5," places.")') inconsistent
      call GOLD_error(WARNING, mesg)
    endif
  endif

! and remap temperature and salinity to layers


  tv%T(is:ie,js:je,:) = tracer_z_init(temp_z(is:ie,js:je,:),-1.0*z_edges_in,zi,nkml,nkbl,missing_value,G%hmask(is:ie,js:je),nz,nlevs)
  tv%S(is:ie,js:je,:) = tracer_z_init(salt_z(is:ie,js:je,:),-1.0*z_edges_in,zi,nkml,nkbl,missing_value,G%hmask(is:ie,js:je),nz,nlevs)

! In case of a problem , use this.

  if (debug_point) then
    do j=js,je ; do i=is,ie
      if (abs(G%geolonh(i,j)-debug_lon).lt.0.25 .and. abs(G%geolath(i,j)-debug_lat) .lt. 0.25) then
        do k=1,kd
          print *,'i,j,klev,T,S,rho=',i,j,temp_z(i,j,k),salt_z(i,j,k),rho_z(i,j,k)
        enddo
        do k=1,nz
          print *,'i,j,klay,T,S,z=',i,j,k,tv%T(i,j,k),tv%S(i,j,k),zi(i,j,k)
        enddo
      endif
    enddo ; enddo
  endif

! Fill land values

  do k=1,nz ; do j=js,je ; do i=is,ie
    if (tv%T(i,j,k).eq.missing_value) then
      tv%T(i,j,k)=temp_land_fill
      tv%S(i,j,k)=salt_land_fill
    endif
  enddo ; enddo ; enddo

! finally adjust to target density

  ks=nkml+nkbl+1

  if (adjust_temperature) then
    call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
            G%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, eos)
  endif

  deallocate(lon_in,lat_in,lon_out,lat_out,x_in,y_in,x_out,y_out)
  deallocate(z_in,z_edges_in)
  deallocate(temp_z,salt_z,rho_z,mask_z)
  deallocate(tmp_in,temp_in,salt_in,mask_in,last_row)
  deallocate(temp_out,salt_out,mask_out,temp_prev,salt_prev)
  deallocate(rho_out,fill,hmask,good,Depth)

  return

end subroutine GOLD_temp_salt_initialize_from_Z

end module GOLD_initialization
